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CHAPTER  1  INTRODUCTION 

1.1  What  is  a  Porous  Medium? 

Almost  anything  can  be  considered  as  porous  depending  on  the  observation  scale 
(Cushman  1997).  In  our  daily  life,  we  usually  think  of  porous  media  as  the  hetero¬ 
geneous  systems  consisting  of  a  solid  matrix  with  fluid-filled  voids.  Dullien  (1992) 
proposed  two  criteria  that  a  material  must  satisfy  in  order  to  qualify  as  a  porous 
medium: 

1.  It  must  contain  relatively  small  spaces,  so-called  pores  or  voids,  free  of  solids, 
embedded  in  the  solid  or  semi-solid  matrix.  The  pores  usually  contain  a  fluid, 
such  as  air,  water,  oil,  etc.,  or  a  mixture  of  different  fluids; 

2.  It  must  be  permeable  to  a  variety  of  fluids,  i.e.,  fluids  should  be  able  to  penetrate 
through  one  face  of  a  septum  made  of  the  material  and  emerge  on  the  other 
side.  In  this  case  one  refers  to  a  "permeable  porous  material  . 

Porous  media  are  ubiquitous  throughout  nature.  Artificial  (man-made)  porous 
media  include  fabric,  paper,  concrete,  cement,  and  brick.  Vegetal  and  animal  bi¬ 
ological  media,  such  as  arteries  and  lungs,  must  convey  fluids  to  transport  oxygen, 
nutrients,  and  wastes.  Geological  porous  media  are  of  interest  in  this  report  because  of 
their  practical  importance  in  oil  recovery,  groundwater  flow,  and  contaminant  trans¬ 
port  in  the  subsurface. 

1.2  Problem  Statement 

Analysis  of  flow  and  transport  processes  through  porous  media  has  many  applica¬ 
tions  in  engineering  systems,  science,  and  technology.  Such  applications  include  the 
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problem  of  seepage,  intrusion  of  sea  water  in  coastal  areas,  contaminant  transport 
through  the  subsurface,  enhanced  oil  recover}',  lubrication,  filtration  of  suspended 
solids,  geothermal  energy  management,  and  nutrition  supply  in  soils.  These  appli¬ 
cations  are  typically  interdisciplinary  and  involve  the  areas  of  civil,  environmental, 
petroleum,  chemical,  mechanical,  and  agricultural  engineering:  geology  and  ground- 
water  hydrology:  food,  soil  and  biochemical  sciences. 

With  the  increasing  sense  of  environmental  awareness,  people  have  interest  in 
the  mechanisms  responsible  for  contaminant  transport  through  groundwater  systems. 
Much  of  the  interest  has  resulted  from  the  enactment  of  new  legislation  which  ad¬ 
dress  the  disposal  of  solid  and  hazardous  wastes  (RCRA)  and  the  cleanup  of  pre¬ 
viously  contaminated  sites  (CERCLA).  The  EPA  Superfund  Innovative  Technology 
Evaluation  (SITE)  Program  is  working  to  expedite  environmental  remediation  by  de¬ 
veloping  promising  new  technologies  for  the  cost-effective  cleanup  of  contaminated 
soil  and  groundwater  (Nyer  1992).  Successful  implementation  of  these  technologies — 
such  as  in  situ  bioremediation,  air  sparging,  soil  vapor  extraction,  surfactant  flushing: 
and  electrokinetic  remediation — requires  a  clear  understanding  of  subsurface  contam¬ 
inant  transport  processes.  To  prevent  the  deterioration  of  groundwater  quality  and 
evaluate  remediation  technologies,  considerable  efforts  have  been  made  toward  the 
development  of  mass  transport  models  for  monitoring,  analyzing,  and  predicting  the 
transport  of  contaminants  through  the  subsurface. 

The  ability  to  model  transport  through  porous  media  has,  however,  been  lim¬ 
ited  by  an  insufficient  understanding  of  physical  processes  at  the  particulate  level. 
Although  the  geometry  of  voids  and  solids  in  porous  media  is  inherently  discrete 
at  the  microscale,  transport  through  porous  media  has  been  traditionally  character¬ 
ized  using  a  macroscopic  continuum  approach  based  on  Darcy’s  law  and  Fick’s  law 
(Ogata  1970:  Bear  1972;  Freeze  and  Cherry  1979).  Numerous  theoretical  and  exper¬ 
imental  studies  have  attempted  to  predict  macroscopic  transport  from  known  media 
properties,  but  knowledge  of  flow  paths  and  mixing  processes  that  occur  within  the 
individual  pores  remains  unsatisfactory.  To  improve  available  continuum  models,  a 
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better  understanding  of  the  fundamental  physics  which  govern  the  flow  and  transport 
processes  at  the  pore-scale  is  required.  One  means  by  which  this  understanding  can 
be  achieved  is  through  the  development  of  pore-scale  models  that  more  accurately 
represent  the  geometry  of  fluid  and  solid  phases  within  porous  media.  Development 
of  such  models  permits  detailed  study  of  flow  and  transport  phenomena  in  assem¬ 
blages  of  particles  with  varying  shape,  angularity,  orientation,  void  ratio,  and  size 
distribution.  This  could  lead  to  important  discoveries  regarding  the  nature  and  in¬ 
terrelationships  of  permeability  and  hydrodynamic  dispersion  tensors  used  in  current 
mass  transport  models. 

In  addition,  better  knowledge  of  pore-scale  flow  and  transport  processes  could 
significantly  improve  our  ability  to  predict  in  situ  liquid-liquid  chemical  and  biologi¬ 
cal  reaction  rates  (surfactant  flushing,  bioremediation),  and  liquid-solid  partitioning 
(adsorption,  retardation).  Such  knowledge  would  aid  our  understanding  of  the  migra¬ 
tion  behavior  of  light  and  dense  non-aqueous  phase  liquids  (LNAPLs  and  DNAPLs) 
through  the  subsurface,  the  mobility  of  which  are  largely  governed  by  surface  tension 
and  pore  geometry  considerations. 

1.3  Objectives  of  Research 

The  principal  objective  of  this  research  is  to  simulate  microscopic  (pore-scale) 
single-phase  flow  and  tracer  transport  processes  through  porous  media,  and  thereby 
gain  insight  with  regard  to  macroscopic  flow  and  transport  behavior.  Rather  than 
employ  a  continuum  approach.  Smoothed  Particle  Hydrodynamics  (SPH)  will  be 
used  to  simulate  fluid  flow  and  solute  transport  processes  within  a  periodic  pore 
network.  SPH  is-well  suited  for  this  study  because  it  is  a  fully  Lagrangian  technique 
in  which  the  numerical  solution  is  achieved  without  a  grid.  Using  this  approach, 
fluid  velocity,  pressure,  contaminant  distributions,  and  volumetric  flow  rates  can  be 
computed.  In  addition,  the  flow  paths  of  individual  fluid  masses  can  be  monitored 
as  they  travel  through  the  void  system  of  the  medium.  Detailed  information  about 
transport  processes  obtained  by  this  approach  would  be  difficult  or  impossible  to 
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observe  experimentally  or  with  many  other  numerical  techniques.  Another  advantage 
of  SPH  is  the  relative  ease  with  which  new  physics  may  be  incorporated  into  the 
formulation.  For  example,  surface  tension  can  be  implemented  into  the  model  to 
simulate  multiphase  transport  (Morris  1999). 

1.4  Outline  of  Report 

Chapter  2  of  this  report  reviews  literature  related  to  the  current  research.  Chap¬ 
ters  3  to  5  present  the  development,  verification,  and  application  of  the  numerical 
model  for  pore-scale  flow,  diffusion,  tracer  convection  and  hydrodynamic  dispersion, 
respectively.  Chapter  6  presents  the  conclusions  and  expected  future  research  using 
the  numerical  model. 


CHAPTER  2  LITERATURE  REVIEW 


Although  practical  hydraulics  had  its  origins  in  antiquity,  scientific  attention  to 
transport  in  particulate  media  began  about  one  hundred  and  fifty  years  ago.  An 
enormous  and  still  rapidly  growing  literature,  including  thousands  ot  articles  and 
manv  monographs,  has  since  been  devoted  to  the  topic  of  flow  and  transport  through 
porous  media.  It  is  not  intended  herein  to  provide  a  comprehensive  review  of  these 
works,  rather,  it  is  the  purpose  of  this  chapter  to  address  the  literature  specifically 
related  to  the  current  research. 

2.1  Models  of  Porous  Media 

A  model  represents  a  useful  simplification  of  complex  reality.  As  the  detailed 
structure  of  natural  porous  media  is  too  complicated  to  describe  mathematically, 
for  years  these  materials  have  been  represented  using  simplified  hypothetical  models 
which  can  be  analytically  treated.  Actually,  any  modeling  of  flow  and  transport 
phenomena  in  a  porous  medium  has  to  include  a  realistic  model  of  the  medium  itself. 
There  have  been  efforts  to  describe  pore-scale  transport  through  spatially  periodic 
porous  media,  network  models,  .fractal  porous  media,  and  reconstructed  porous  media. 
The  following  sections  outline  the  progress  which  has  been  made  in  these  areas. 

2.1.1  Spatially  Periodic  Model 

A  spatially  periodic  model  views  the  porous  medium  as  having  a  spatially  periodic 
structure.  This  is  based  on  translational  symmetry,  which  is  a  classical  concept 
that  many  materials  look  much  the  same  at.  different  locations.  Spatially  periodic 
structure  is  the  simplest  structure  that  can  be  seen  in  nature  and  was  historically  the 
first  to  be  studied.  Spatially  periodic  media  are  generated  by  infinite  sequences  of 
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repeating  unit  cells  in  one,  two.  or  three  dimensions.  For  a  three-dimensional  problem, 
periodic  media  can  be  constructed  following  the  work  of  Adler  (1992).  Consider  a 


Figure  2.1:  Unit  cell  of  a  spatially  periodic  structure. 

parallelepiped  unit  cell  whose  sides  are  characterized  by  three  vectors  Ii,  I2,  and  I3 
(Figure  2.1).  The  content  of  this  cell  is  arbitrary.  It  can  be  composed  of  particles,  of 
a  network  of  capillary  tubes,  or  of  any  other  solid  or  void  organization.  Its  volume 
I  unit  ceil  is  assumed  to  be  strictly  positive, 

I  unit  cell  —  1^1  '  ^2  ^3 1  ?  (2*1) 

which  is  equivalent  to  assuming  that  these  three  vectors  are  linearly  independent.  An 
infinite  spatially  periodic  medium  is  obtained  by  translating  this  unit  cell  by  all  linear 
combinations  R„  of  the  three  vectors  lx,  I2,  and  I3, 

Rn  =  nilj  +  n2 12  +  n-3l3i  (2-2) 

where  the  trio  (rij,  n2.  n3)  are  integers.  By  doing  this,  infinite  unit  cells  are  simply 
juxtaposed  to  one  another  to  form  a  spatially  periodic  medium.  If  R„  is  defined  as  the 


centroid  of  each  unit  cell,  the  global  position  vector  R  anywhere  within  the  infinite 
medium  mav  be  identified  as. 


R  =  Rn  +  r.  (2.3) 

where  r  is  the  local  position  vector  originating  at  the  centroid  of  each  unit  cell  (Figure 
2.1).  Figure  2.2  depicts  an  irregular  but  spatially  periodic  porous  medium  in  two 
dimensions. 


Figure  2.2:  A  spatially  periodic  porous  medium. 


A  spatially  periodic  medium  is  characterized  by  three  physical  length  scales.  The 
first  is  the  microscopic  length  scale  Cm  of  the  characteristic  size  of  the  particles  over 
which  the  local  interstitial  fields  (e.g.,  velocity  field)  vary  due  to  the  local  boundary 
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conditions  satisfied  on  the  particle  surfaces.  The  second  is  the  Darcy  scale  of 
the  order  of  the  size  of  the  unit  cell  over  which  the  mean  or  average  fields  vary 
sensibly.  People  are  usually  interested  in  physically  describing  the  gross  or  average 
transport  process  occurring  at  the  Darcy  scale.  Quantities  defined  at  the  Darcy  scale 
represent,  in  some  sense,  averages  of  comparable  quantities  defined  at  each  point  of 
the  interstices.  It  is  necessary  to  require. 

CD^Cm:  (2.4) 

such  that  a  volume  element  of  0(C3D)  contains  a  representative  number  of  particles 
and  encapsulates  enough  degree  of  heterogeneity  to  render  the  various  meaningful 
averages  formed  from  integration  of  the  local  interstitial  fields.  The  last  is  the  macro¬ 
scopic  length  scale  C  which  corresponds  to  a  characteristic  linear  dimension  of  the 
external  boundaries  of  the  porous  medium  under  consideration.  For  example,  the  di¬ 
mension  of  a  specimen  of  porous  media  upon  which  an  experiment  is  being  performed. 
In  order  to  model  a  porous  medium  as  spatially  periodic,  it  is  assumed  that. 

C  =  oo  ■  and  L  )§>  Co,  (2-5) 

so  that  some  average  fields  defined  on  the  Darcy  scale  remain  sensibly  constant, 
such  as  pressure  gradient  and  the  resulting  average  discharge  velocity.  Similarly,  the 
average  concentration  gradient  and  the  corresponding  mass  flux  are  sensibly  constant 
over  Cd- 

Some  theoretical  and  numerical  research  has  been  devoted  to  transport  through 
spatially  periodic  porous  media.  By  assigning  a  given  solid  configuration  to  the  unit 
cell,  it  is  possible  to  study  the  problem  at  the  pore-scale.  As  this  report  also  assumes 
a  porous  medium  has  spatial  periodicity,  the  related  literature  will  be  reviewed  in 
detail  in  section  2.4. 

2.1.2  Network  Model 

The  use  of  networks  of  capillary  tubes  to  model  porous  media  was  first  suggested 
by  Fatt  (1956).  The  work  was  based  on  the  idea  that  pore  space  may  be  represented 
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as  an  interconnected  network  of  capillary  tubes  whose  radii  represent  the  dimensions 
of  the  pores  within  a  porous  medium.  The  premise  of  the  network  model  is  that  the 
void  space  of  a  porous  medium  can  be  represented  by  a  graph  of  connected  nodes. 
The  nodes  in  the  graph  correspond  to  pore  bodies  and  the  links  that  connect  the 
nodes  correspond  to  pore  throats.  Such  a  graph  preserves  the  essential  topology  of 
the  pore  spaces.  In  addition,  sizes  can  be  assigned  to  the  nodes  and  links  according 
to  a  chosen  size  distribution  of  pore  bodies  and  throats  (Bryant  et  al.  1993). 

The  pore-scale  network  approach  for  investigating  the  nature  of  fluid  flow  has 
been  developed  extensively  in  the  petroleum  engineering  (Chatzis  and  Dullien  1977: 
Larson  et  al.  1981),  and  also  recently  in  the  fields  of  hydrology  and  soil  physics 
(Ferrand  and  Celia  1992;  Reeves  and  Celia  1996;  Rajaram  et  al.  1997).  Although 
network  models  can  replicate  highly  disordered  geometry  of  fluid  phases  in  principal, 
in  practice,  they  have  been  constructed  using  various  assumptions  concerning  pore 
structure.  Mathematical  network  models  are  used  in  computer  simulations  while 
physical  network  models  have  been  developed  for  flow  visualization  studies  (Wan 
et  al.  1996). 

2.1.3  Fractal  Model 

Mandelbrot  (1982)  coined  the  name  “fractal”  from  the  Latin  fractus  which  de¬ 
scribes  the  appearance  of  a  broken  stone  and  popularized  the  concept  of  fractals.  It 
has  been  found  in  many  fields  that  fractals  could  mimic  the  geometry  of  the  real 
world  better  than  classic  geometry.  Fractals  provide  a  powerful  concept  for  physicists 
to  solve  problems  characterized  by  the  simultaneous  existence  of  very  different  scales. 
In  recent  years,  many  complex  behaviors  of  a  wide  variety  of  phenomena  have  been 
quantitatively  characterized  using  the  concept  of  fractals. 

A  number  of  authors  have  studied  the  fractal  character  of  porous  media  (Avnir 
et  al.  1985;  W’ong  et  al.  1986;  Jacquin  and  Adler  1987;  Thompson  et  al.  1987).  Frac¬ 
tal  porous  media  are  highly  self-organized  in  that  many  static  and  dynamic  properties 
have  infinite  correlation  lengths.  The  problem  scale  of  fractal  porous  media  evolves 
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continuously  from  pore  to  reservoir.  Constitutive  variables  in  such  a  system  are  wave 
vector  and  frequency  dependent,  and  result  from  non-local  theories  (Cushman  1990: 
Cushman  and  Ginn  1993:  Cushman  et  al.  1994:  Cushman  1997). 

2.1.4  Reconstructed  Porous  Media 

Reconstruction  of  actual  porous  media  has  been  done  using  computer  microto¬ 
mography  (CMT)  (Spanne  et  al.  1994:  Peyton  et  al.  1994:  Schwartz  et  al.  1994)  and 
digital  image  processing  (Koplik  et  al.  1984).  Spanne  et  al.  (1994)  applied  computer 
microtomography  to  sandstone  to  determine  the  geometrical  structure  of  the  port's 
to  a  resolution  of  10  pm.  Reconstructed  porous  media  may  mimic  more  closely  the 
geometry  of  real  media  and  generate  numerical  samples  with  desired  properties. 

Adler  (1992)  presented  three  steps  for  the  study  of  transport  in  reconstructed 
porous  media.  The  first  involves  the  measurement  of  any  salient  geometric  features. 
The  second  step  is  the  reconstruction  process.  Samples  of  porous  media  are  generated 
in  such  a  way  that,  on  average,  they  possess  the  same  statistical  properties  as  the  real 
samples  that,  they  are  intended  to  mimic.  For  the  last  step,  all  transport  phenomena 
are  studied  analytically  or  numerically.  Reconstruction  of  porous  media  can  be  done 
on  pore-  or  field-scale,  depending  on  the  problem. 

2.2  Single-Phase  Flow  Through  Porous  Media 

2.2.1  Introduction 

Henry  P.  G.  Darcy  published  studies  on  the  development  of  the  water  supply 
systems  of  Dijon  in  1856.  The  law  he  discovered,  namely,  that  the  rate  of  flow  is 
proportional  to  the  total  head  drop  through  a  bed  of  fine  particles,  bears  his  name 
and  is  widely  employed  for  investigating  the  behavior  of  all  types  of  fluid  flow  through 
porous  media. 

For  isotropic  media,  Darcy’s  law  is  written  in  differential  form  as, 


v  =  -KVP  =  -I<(Vp  +  pg) 


(2.6) 
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where  v  denotes  the  filter  velocity  (i.e.,  Darcy  velocity,  discharge  velocity,  or  specific 
discharge)  and  VP  is  the  pressure  gradient  vector.  I\  is  the  permeability,  p  is  the 
hydrostatic  pressure,  p  is  the  fluid  density,  g  is  the  gravitational  acceleration,  and  the 
total  pressure  P  is  defined  as, 

P  =  p  +  pgz ,  (2.7) 

where  2  is  the  distance  measured  vertically  upward  from  an  arbitrarily  chosen  datum 
level,  and  g  is  the  gravitational  constant.  Usually,  P  is  measured  by  a  piezometer 
and  is  indicated  as  the  "piezometric  head"  or  “total  head  h. 

h=  — =  2- +  z,  (2.S) 

pg  pg 

T) 

which  is  the  sum  of  the  elevation  head  z  and  the  pressure  head  ^ . 

The  permeability  K  depends  on  both  the  medium  and  the  fluid.  Nutting  (1930) 
separated  the  influence  of  the  porous  medium  from  that  of  the  fluid  and  stated  that. 

K  =  -,  (2-9) 

P 

where  //  is  the  dynamic  viscosity  of  the  fluid  and  k  is  called  the  specific,  absolute,  or 
intrinsic  permeability  of  the  medium  which  has  dimensions  of  length  squared.  This 
concept  was  popularized  by  Wyckoff  et  al.  (1933).  For  simplicity,  k  will  be  called 
permeability  in  this  report. 

Flow  in  porous  media  takes  place  through  flow  channels,  each  having  a  local 
distribution  of  pore  velocity.  The  average  pore  velocity,  on  the  whole,  must  be  larger 
than  the  Darcy  velocity  owing  to  the  reduced  space  available  for  flow,  as  compared 
with  the  bulk  volume  of  the  porous  medium  on  which  the  Darcy  velocity  is  defined. 
A  commonly  accepted  hypothesis  for  the  connection  between  pore  velocity  and  Darcy 
velocity  is  the  Dupuit-Forchheimer  assumption  (Scheidegger  1974), 

v  =  vsn,  (2.10) 

where  v5  is  the  pore  or  seepage  velocity  and  n  is  the  porosity  of  the  medium.  It 
should  be  noted  that  the  seepage  velocity,  like  Darcy  velocity,  is  not  a  true  physical 
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velocity  because  the  actual  velocity  of  the  fluid  must  be  expected  to  fluctuate  grossly 
within  one  flow  channel  and  from  one  flow  channel  to  another,  and  because  of  the 
tortuosity  of  the  channels. 

In  groundwater  hydrology  and  soil  mechanics,  the  fluid  of  interest  is  usually  water. 
Therefore,  a  "hydraulic  conductivity"  kn  is  defined  as. 

*„  =  h!£.  (2.m 

and  Darcy's  law  is  written  as, 

v  =  -kHVh=  -kH  i,  (2.12) 

where  i  =  V/i  is  the  hydraulic  gradient. 

The  concept  of  permeability  permits  a  phenomenological  description  of  flow  through 
porous  media.  However,  an  actual  understanding  of  the  phenomena  can  be  obtained 
only  if  the  concept  of  permeability  can  be  reduced  to  more  fundamental  physical  prin¬ 
ciples.  It  is  intuitively  clear  that  permeability  is  linked  with  other  properties  of  porous 
media.  Many  attempts  have  been  made  to  establish  correlations  between  permeabil¬ 
ity  and  other  properties  of  porous  media  with  the  help  of  various  models.  These 
approaches  may  be  categorized  in  several  different  ways.  Here,  two  fundamentally 
different  approaches  are  distinguished:  in  one,  the  flow  inside  conduits  is  analyzed; 
in  the  other,  the  flow  around  solid  objects  immersed  in  the  fluid  is  considered. 

2.2.2  Conduit  Flow  Model 

Conduit  flow  does  not  account  for  the  fact  that  different  pores  are  intercon¬ 
nected  with  each  other,  so  conduit  flow  models  are  inherently  one-dimensional  models. 
Among  them,  tRe  Kozeny-Carman  model  is  generally  more  popular  than  the  rest. 

2. 2. 2.1  Kozeny-Carman  Model 

The  Kozeny-Carman  approach  is  also  called  the  “hydraulic  radius  theory”  and 
was  developed  for  creeping  flow.  In  the  Kozeny-Carman  theory  (Kozeny  1927;  Car¬ 
man  1937:  Carman  1938a;  Carman  1938b;  Carman  1939;  Carman  1956),  the  porous 
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medium  is  assumed  to  be  equivalent  to  a  conduit,  the  cross  section  of  which  mac  hace 
a  complicated  shape  but,  on  an  average,  a  constant  area. 

A  Hagen-Poiseuille  type  equation  is  assumed  to  give  the  average  seepage  velocity 
in  the  flow  channels, 

„=.^_AL  r2.n1 

L,  C,  16/i' 

where  Le  is  the  effective  path  length  of  flow.  D  is  the  flow  channel  diameter,  and  C, 
is  a  shape  factor. 

In  addition  to  the  Dupuit-Forchheimer  assumption  in  Equation  2.10.  Carman 
related  Darcy  velocity  with  the  seepage  velocity  as. 


v 


S 


V_  Le 
nl’ 


(2.14) 


which  corrects  for  the  fact  that  a  hypothetical  fluid  particle  used  in  macroscopic  flow- 


equations  and  flow-ing  with  velocity  v  covers  a  path  length  L  in  the  same  time  as  an 
actual  fluid  particle  flow-ing  w-ith  velocity  vs  covers  an  average  effective  path  length 
Le,  as  show-n  in  Figure  2.3. 


Combination  of  Equations  2.13  and  2.14  gives, 
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where  D  is  assumed  to  be  four  times  the  hydraulic  radius, 

void  volume  of  medium 


D 


4  n 


>.16) 


surface  area  of  channels  in  medium  S0(l  —  n) 
and  So  is  the  specific  surface  area  based  on  the  solid  volume.  By  combining  Equations 
2.15  and  2.16,  the  usual  form  of  the  Kozeny-Carman  equation  for  permeability  is. 


k  = 


n 


(2.17) 


C 's  ("Tf )  _  n)2So 

where  ^  is  the  tortuosity  factor  T.  A  combined  factor  k'  =  CST^  =  Cs  K  t ho 

Kozeny-Carman  constant.  The  Kozeny-Carman  equation  indicates  the  overwhelming 
importance  of  pore  size  (through  the  So  term)  in  determining  the  permeability  of  a 
porous  medium. 

According  to  Carman,  Cs  lies  within  the  range  2.0  to  3.0  with  a  likely  average 
value  of  2.5.  He  also  suggested  that  T2  has  a  value  of  about  2.0  in  all  unconsolidated 
porous  aggregates  like  packed  beds,  which  results  in  k'  =  5.0.  Much  evidence  that  is 
now  available  certainly  suggests  that  for  unoriented  particle  aggregates  in  the  porosity 
range  of  0.35  to  0.70,  k'  =  5.0  ±  10%  according  to  Wyllie  and  Spangler  (1952)  and 
Wvllie  and  Gregory  (1955).  They  also  found  k'  is  dependent  on  the  porosity,  particle 
shape,  and  particle  orientation.  For  aggregates  of  fibers  at  very  large  porosities  (n  > 
0.84),  k'  was  found  to  increase  rapidly  with  the  porosity. 

Although  Scheidegger  (1974)  has  strong  criticism  of  the  Kozeny-Carman  equa¬ 
tion,  some  scholars  consider  the  equation  to  be  approximately  valid  for  certain  soil 
types  (Dullien  1992;  Mitchell  1993).  The  Kozeny-Carman  equation  can  be  applied  to 
uniformly  graded  sands  and  some  silts,  although  serious  discrepancies  are  found  when 
the  equation  is  applied  to  clays.  The  main  reason  of  the  discrepancy  is  that  clays  do 
not  contain  approximately  uniform  pore  sizes.  Particles  in  clays  are  grouped  in  clus¬ 
ters  or  aggregates  that  produce  large  intercluster  pores  and  small  intracluster  pores. 
It  is  implicit  in  the  derivation  of  the  Kozeny-Carman  equation  that  there  shall  be  no 
large  and  small  pores  in  parallel  contributing  to  flow.  If  such  a  condition  exists,  the 
contribution  to  flow  made  by  the  larger  pores  is  disproportional  to  the  corresponding 
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effect  on  specific  surface  area.  The  latter  quantity  is  dominated  by  the  dimensions  of 
the  small  pores. 

2. 2. 2. 2  Capillary  Model 

The  capillary  permeability  model  is  the  simplest  approach  based  on  the  idea  of 
conduit  flow.  Scheidegger  (1974)  distinguished  among  straight,  parallel,  and  serial 
type  capillary  models. 

The  straight  capillary  model  (Figure  2.4)  represents  a  porous  medium  as  a  bundle 
of  identical  and  parallel  straight  capillaries  of  uniform  diameter  D.  The  total  volume 


Figure  2.4:  Straight  parallel  capillary  model. 


flow  rate  Q  through  a  capillary  is  given  by  the  law  of  Hagen-Poiseuille, 


7rP4  dP_ 
12  Sfj,  dx 


(2.18) 


If  there  are  m  such  capillaries  per  unit  area  of  cross  section  of  the  model,  the  Darcy 


velocity  will  be, 


and, 


mirD 4  dP 
128 //  dx  ’ 

rmiD4 


(2.19) 


(2.20) 
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As  the  model  has  a  porosity  of. 


1  r-\2 

n = -m~D  . 
4 


(2.21) 


the  permeability  is. 


k  = 


nD2 

HT' 


12.22) 


which  is  a  variant  of  Kozenv-Carman  equation.  Equation  2.22  also  bears  similarity 
to  the  empirical  Hazen  equation. 


k  —  C  D 2 


2.23) 


where  Cm  is  an  empirical  parameter  and  D10  is  the  effective  grain  size  of  a  loose  clean 
sand. 

In  the  straight  capillary  model,  there  is  no  flow  orthogonal  to  the  capillaries. 
A  general  parallel  model  of  capillaries  with  different  diameters  has  one-third  of  the 
capillaries  in  each  of  three  spatial  dimensions  and  arrives  at  the  permeability  as. 


=-  r 

96  Jo 


D2a(D)dD. 


(2.24) 


The  mean  square  diameter  is  evaluated  according  to  a  pore  size  distribution  function 
a[D)  and  96  is  substituted  for  32  in  Equation  2.22  because  only  ^  of  the  capillaries 
are  oriented  in  each  direction.  Equation  2.24  is  sensitive  to  error  at  the  extrema  of  the 
pore  size  distribution  function  corresponding  to  the  largest  pore  sizes  (Scheidegger 
1974). 

In  the  serial  type  capillary  model  (Figure  2.5),  the  pore  network  is  approximated 
by  three  identical  sets  of  tortuous  channels.  Each  channel  is  assumed  to  consist  of 
segments  having  different  diameters  distributed  according  to  a  pore  size  distribution 
function  a(D).  The  permeability  of  this  model  is, 


k  = 


n 


(/ 


a(fl) 

D2 


dD 


96 


/ 


<*{D) 

D6 


dD 


(2.25) 


The  serial  type  model  tends  to  underestimate  permeabilities  and  is  very  sensitive  to 
uncertainties  at  the  extrema  of  the  pore  size  distribution  function  corresponding  to 
the  smallest  pores  (Dullien  1992). 


2.2.3  Drag  Flow  Model 

As  discussed  by  Scheidegger  (1974),  Emersleben  (1925)  first  proposed  the  'drag 
theory  of  permeability”,  which  was  different  from  the  theory  of  Kozeny.  This  approach 
considers  flow  around  submerged  objects.  The  walls  of  the  pores  are  modeled  as 
obst  acles  in  an  otherwise  straight  flow  of  viscous  fluid.  Fluid  drag  forces  are  estimated 
from  the  Navier-Stokes  equations,  and  the  sum  of  drag  forces  is  assumed  to  give  the 
resistance  of  porous  medium  to  flow.  The  subject  of  low  Reynolds  number  flow  around 
submerged  objects  is  discussed  in  detail  by  Happel  and  Brenner  (1965). 

Iberall  (1950)  specified  the  obstacles  as  cylinders.  He  considered  the  permeability 
of  a  random  distribution  of  circular  cylinders  and  showed  for  low  Reynolds  number 
Re  that, 

3  nD 2  2  —  In  Re 
16  1  —  n  4  —  In  Re ' 


(2.26) 
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where, 


pvsD 


(2.27) 


D  is  the  diameter  of  the  cylinders,  and  vs  is  t he  average  pore  velocity  ( i.e. .  seepage 
velocity).  It  is  noted  from  Equation  2.26  that  permeability  is  not  constant  but  varies 
with  the  Reynolds  number,  i.e.,  the  flow  velocity,  although  the  derivation  was  actu¬ 
ally  based  on  the  assumption  of  low  Reynolds  number  flow.  This  slow  variation  of 
permeability  with  flow  is  quite  characteristic  of  many  instances  in  which  the  How  is 
nominally  viscous. 

Brinkman  (1947,  1948)  assumed  the  obstacles  in  the  fluid  were  spheres  of  diameter 
D  held  in  position  by  external  forces.  He  extended  the  Stokes  drag  force  on  a  sphere  (a 
sphere  placed  in  an  infinite  flow  domain)  to  include  the  effect  of  neighboring  spheres 
and  superimposed  Stokes  and  Darcy  flow  to  obtain. 


VP  =  -yV  +  gV2v. 

K 

The  relationship  derived  by  Brinkman  for  the  permeability  is, 


(2.2S) 


,  D 2  (  4 

k  —  —  3  ~b - 

72  \  1  —  n 


1  —  n 


(2.29) 


Tnis  equation  gives  k  =  0  for  n  =  which  makes  it  unsuitable  for  relatively 
low  porosities  (Happel  and  Brenner  1965).  Lundgren  (1972)  proposed  a  generalized 


Brinkman  equation, 


72  F{n) 


4  /  8 

3  +  l - 3\/l - 

1  —  n  VI  —  n 


Equation  2.30  includes  an  “effective  viscosity”  T(n)  defined  as, 

.  .  47:  a2R 2 

3  (1  —  n)F(a2R2,aR), 


(2.30) 


(2.31) 


33  +  vni"3 

4  13' 


(2.32) 
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r  =  D  and  F[a2R2,aR)  is  a  complicated  expression  involving  Bessel  functions 
and  Legendre  polynomials.  However,  the  modified  Brinkman  equation  behaves  much 
differently  from  other  predictions  for  n  <  0.7  (Dullien  1992). 

The  resistance  of  arbitrarily  shaped  objects  in  a  steady  state  stream  of  a  viscous 
fluid  can  be  calculated  by  solving  the  Navier-Stokes  equations  for  the  appropriate 
boundary  conditions.  However,  analytical  solutions  for  such  problems  are  very  lim¬ 
ited. 

2.2.4  Anisotropic  Porous  Media  and  Permeability  Tensor 

Most  porous  media  are  anisotropic.  Anisotropic  permeability  can  result  from  both 
preferred  orientation  of  elongated  or  platy  particles  and  stratification  of  soil  deposits 
(Mitchell  1993).  While  the  pressure  gradient  and  Darcy  velocity  vectors  are  parallel 
in  isotropic  media,  their  directions  are  generally  non-parallel  in  anisotropic  media. 
Consequently,  the  value  of  permeability  depends  on  the  direction  of  measurement  in 
an  anisotropic  porous  medium. 

In  anisotropic  media,  Darcy’s  law  is  written  as, 

v  =  — — VP,  (2.33) 

where  the  permeability  k  takes  the  form  of  a  second-order  tensor  in  a  x-y-z  Cartesian 
coordinate  system, 


Considerable  effort  has  been  expended  to  prove  the  k  tensor  is  symmetric,  i.e.,  kxy  = 
kyx,  etc.  (Szabo  1968;  Whitaker  1969;  Gum  et  al.  1971;  Case  and  Cochran  1972). 
The  physical  meaning  of  these  equalities  is  that  the  permeability  values  in  opposite 
directions  in  the  medium  are  the  same,  so  that  the  k  tensor  changes  the  axes  without 
deformation  of  the  system.  In  exceptional  cases,  where  the  pores  of  a  medium  form  a 
symmetric  helicoidal  arrangement  and  the  velocity  field  is  forced  to  have  rotation,  the 


k  = 


kXx  kXy  kXZ 
k  k 

^ yy  a yz 

kzx  kZy  kzz 


(2.34) 
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tensor  is  not  symmetric  (Liakopoulos  1965).  It  is  generally  assumed  that  anisotropic 
porous  media  are  “orthotropic”,  i.e..  they  have  three  mutually  orthogonal  principal 
axes.  For  an  orthotropic  medium,  rotation  of  the  coordinate  system  will  produce  a 
diagonal  k  matrix  when  the  coordinate  axes  are  aligned  with  the  principal  axes  of 
the  medium.  In  this  case.  VP  and  v  are  parallel  and  Darcy's  law  becomes. 


kx  8P 

(2.35) 

v'  ~~Ji  d r 

k2  dP 

V 2  = - -X-, 

p  o2 

(2.36) 

ki  dP 

n  dz  ’ 

(2.37) 

where  1,  2,  and  3  denote  the  three  principal  directions  of  the  permeability  tensor, 
and  Vi  and  ki  ( i  —  1,2,3)  are  the  Darcy  velocities  and  permeabilities  in  the  three 
principal  directions,  respectively. 

In  practice,  the  permeability  is  usually  measured  along  a  direction  n.  which  may 
not  necessarily  coincide  with  the  principal  permeability  directions.  Darcy’s  law  in 
this  direction  is  written  as, 


kn  dP 

Vn  = - -X- 

/i  an 


(2.38) 


pjP 

where  vn  and  are  the  projections  of  v  and  VP  in  the  n  direction,  respectively, 
i.e.. 


Vn  =  n  •  V, 


(2.39) 


dP 

—  =  n  •  VP,  (2.40) 

on 

and  kn  is  called  the  directional  permeability.  The  directional  permeability  may  not 
necessarily  be  equal  to  the  projection  of  k  in  the  n  direction  (Case  1971). 

Scheidegger  (1974)  investigated  the  relation  between  directional  permeability  and 
permeability  tensor.  He  distinguished  two  possible  cases.  Case  1:  the  velocity  is 
measured  directly  and  the  component  of  the  pressure  gradient  in  the  direction  of  the 
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velocity  is  used  in  Darcy’s  law;  Case  2:  the  velocity  in  the  direction  of  the  pressure 
gradient  is  measured  and  used  in  Darcy's  law. 

In  the  first  case,  from  Equations  2. 38.  2.39  and  2.40. 

■  “TO 

Operating  on  the  tensor  form  of  Darcy’s  law  in  Equation  2.33  leads  to. 

-VP  =  p  k_1-v.  (2.42) 

where  k_1  is  the  inverse  of  k.  Combining  Equations  2.41  and  2.42  gives. 

k  =  v  = _ ! _ .  (2.4-3) 

n  n-k-J-v  rik-’n 

If  the  principal  axes  of  the  permeability  tensor  are  chosen  as  coordinate  axes  and  the 
n  direction  makes  angles  a.  / 3 ,  and  7  with  the  principal  axes.  Equation  2.43  becomes. 


1  cos2  q  cos2  p  cos2  ~ 

—  = - 1 - -  h - 

kn  k\  ko  k^ 


(2.44) 


In  two  dimensions,  according  to  Equation  2.44,  the  plot  of  A*„  produces  an  ellipse 
whose  axes  are  in  the  principal  directions  of  permeability  with  lengths  of  axes  equal 
to  twice  the  square  root  of  the  principal  permeability  values  (Figure  2.6). 

In  the  second  case,  however,  pressure  gradient  is  in  the  n  direction,  i.e.. 

VP  =  ^n,  (2.45) 

on 


1  .  ,  „  1  „  ,  dP  „ 

vn  =  n  •  v  = - n  •  kVP  =  — n  •  k— n. 

jj,  p  on 


Comparison  of  Equations  2.46  and  2.38  yields, 

kn  —  n  •  k  •  n. 


(2.46) 


(2.47) 


If  the  principal  cixes  of  the  permeability  tensor  are  coordinate  axes,  a  similar  treatment 
as  that  used  in  Case  1  results  in  following  expression, 

kn  =  ki  cos2  a  +  /c2  cos2  (3  +  k3  cos2  7.  (2.48) 
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Equation  2.48  defines  an  ellipse  if  kn  2  is  plotted  for  various  n  directions  in  two 
dimensions  (Figure  2.7). 

Given  the  first  ellipse  of  permeability,  the  direction  of  pressure  gradient  can  be 
graphically  determined  for  a  given  velocity,  while  the  direction  of  velocity  for  a  gi\en 
pressure  gradient  can  be  determined  from  the  second  ellipse  of  permeability  (Li- 
akopoulos  1965;  Case  and  Cochran  1972).  The  construction  procedure  for  obtaining 
the  direction  of  the  velocity  vector  for  a  known  direction  of  pressure  gradient  is  (Li- 
akopoulos  1965): 

■  1.  Construct  the  permeability  ellipse  in  the  directions  of  its  principal  axes  with 

semi-axes  equal  to  the  inverse  square  root  of  the  principal  permeabilities,  i.e.. 
the  second  ellipse  of  permeability; 

2.  Draw  the  direction  of  the  pressure  gradient  vector  passing  through  the  center 
of  the  ellipse; 

3.  Draw  the  tangent  plane  at  the  point  where  the  pressure  gradient  vector  pene¬ 
trates  the  surface  of  the  ellipse; 

4.  Draw'  the  normal  to  the  tangent  plane  at  the  penetration  point.  This  is  the 
direction  at  which  the  flow'  w'ill  take  place  (Figure  2.7). 

A  similar  procedure  can  be  followred  to  determine  the  direction  of  pressure  gradient 
using  the  first  ellipse  of  permeability  (Figure  2.6).  More  details  can  be  found  in 
Liakopoulos  (1965). 

There  is  a  standard  method  for  the  transformation  of  a  symmetric  tensor  in  any 
arbitrary  Cartesian  coordinate  system  to  one  in  w'hich  the  tensor  is  diagonal.  For  a 
two-dimensional  anisotropic  porous  medium,  if  kxx,  kyy,  and  kxy  are  known  in  the  x-y 
Cartesian  coordinate  system  (Figure  2.8),  the  values  of  the  principal  permeabilities 
k\  and  k2  and  the  angle  6  (0°  <  9  <  180°)  between  the  direction  of  k\  and  the  positive 
x  direction  can  be  determined  by  (Szabo  1968), 
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Figure  2.8:  x-y  Cartesian  coordinate  system  and  principal  tensor  axes. 

2.2.5  Limitations  of  Darcy’s  Law 

Non-Darcy  flow  behavior  was  reported  as  early  as  1898  (Mitchell  1993).  The 
following  hypotheses  have  been  proposed  to  account  for  non-linearity  between  flow 
velocity  and  gradient  (Scheidegger  1974;  Mitchell  1993):  (1)  non-Newtonian  behavior 
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of  the  fluid  itself:  (2)  particle  migrations  that  cause  blocking  and  unblocking  of  pores: 
(3)  local  consolidation  and  swelling  of  the  soil:  (4)  molecular  or  ionic  effects:  and  (5) 
turbulent  effects  associated  with  high  flow  rates. 

The  linearity  of  Darcy's  law  can  be  derived  from  the  Navier-Stokes  equations  if 
the  inertial  and  unsteady  terms  are  neglected.  Ignoring  the  inertial  term  restricts  the 
analysis  to  creeping  (laminar)  flow  where  viscous  forces  dominate  in  the  flow  regime. 
The  range  of  validity  of  Darcy's  law  is  expressed  in  terms  of  Re. 


where  £  is  a  characteristic  length  and  V  is  a  characteristic  velocity  of  the  flow  system. 
£  is  generally  the  mean  pore  size,  mean  grain  size,  or  \/k.  Re  is  a  measure  of  the 
ratio  between  the  inertial  and  viscous  forces  of  the  flow.  In  any  event,  it  is  generally 
believed  that  Darcy’s  law  is  valid  as  long  as  Re  <  5. 

Darcv’s  law  can  be  modified  to  account  for  non-linear  flow  behavior  at  high  ve¬ 


locities.  According  to  Scheidegger  (1974),  Forchheimer  (1901)  suggested  that  Darcy  s 
law  include  a  second-order  velocity  term, 


dP 

dx 


=  QiV  +  02l>2, 


(2.52) 


where  aj  and  a2  are  constants.  To  obtain  a  better  representation  of  experimental 
data,  a  cubic  term  was  added  later  (Scheidegger  1974), 


dP 

dx 


=  CL\V  +  0.2V2  +  G3U3. 


(2.53) 


At  low  velocities,  v 2  and  v 3  are  much  less  than  v  and  Equation  2.53  reduces  to 
Darcy’s  law.  Tlie  Forchheimer  Equation  has  been  generalized  further  to  contain  a 
time-dependent  term  (Polubarinova-Kochina  1952), 


dP 

dx 


—  a\v  +  CL2V 


2 


dv 

+a3di- 


(2.54) 


Other  heuristic  corrections  to  Darcy's  law  can  be  found  in  Scheidegger  (1974)  and 
Dullien  (1992). 
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2.3  Hydrodynamic  Dispersion  in  Porous  Media 

2.3.1  Introduction  and  Convection-Dispersion  Equation 

According  to  Bear  (1972),  one  of  the  earliest  observations  of  dispersion  was  re¬ 
ported  by  Slichter  in  1905  who  used  an  electrolyte  as  a  tracer  in  studying  the  move¬ 
ment  of  groundwater.  A  tracer  is  a  chemical  compound  that  does  not  affect  the 
density  or  the  viscosity  of  the  fluid  phase  in  which  it  is  diluted  and  does  not  sorb 
to  the  solid  phase  of  the  porous  medium  during  transport.  Slichter  observed  that 
at  an  observation  well  downstream  of  the  injection  point,  the  tracer  concentration 
increased  gradually  with  time.  It  was  also  observed  that  in  an  uniform  flow  field 
the  tracer  advanced  in  the  direction  of  the  flow  in  a  pear-like  shape  that  became 
longer  and  wider  as  it  advanced.  Slichter  explained  this  phenomena  by  noting  that, 
for  flow  through  capillary  tubes,  fluid  velocity  varies  across  the  cross  section  of  each 
tube  and  that,  because  a  porous  medium  is  composed  of  a  great  number  of  these 
tubes,  the  combined  effect  is  likely  to  cause  the  mixing  that  he  observed.  This  tracer 
spreading  phenomenon  is  called  hydrodynamic  dispersion  in  porous  media.  It  is  a 
nonsteady,  irreversible  process  in  which  the  tracer  mass  mixes  with  progressively 
larger  volumes  of  fluid.  If  the  tracer  initially  occupies  a  separate  region  in  the  flow 
field,  an  ever-widening  transition  zone  is  created  due  to  dispersion,  across  which  the 
tracer  concentration  varies  from  that  of  the  tracer  liquid  to  that  of  the  carrier  liquid. 
In  practice,  the  displacement  of  two  miscible  fluids  may  be  considered  as  tracer  flow 
when  the  fluids  have  the  same  densities  and  viscosities,  and  there  is  no  volume  change 
associated  with  mixing. 

Contaminant  transport  in  the  subsurface  is  generally  viewed  as  the  net  effect  of 
two  processes,  convection  and  dispersion.  Convection  is  transport  by  the  average  mo¬ 
tion  of  the  fluid  and  is  the  primary  mechanism  responsible  for  contaminant  migration 
in  many  aquifers.  Dispersion  is  spreading  of  the  convective  contaminant  front  due  to 
the  movements  of  individual  contaminant  particles  through  the  pores  and  the  various 
physical  and  chemical  phenomena  that  take  place  within  the  pores.  Mechanisms  that 
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contribute  to  dispersion  include  molecular  diffusion,  mechanical  dispersion,  bound 
ary  layer  or  film  diffusion,  and  intraparticle  diffusion.  Currently,  one  of  the  most 
widely  held  tenets  is  the  assumption  that  hydrodynamic  dispersion,  which  results 
from  molecular  diffusion  and  mechanical  dispersion,  constitutes  the  total  dispersion. 


2. 3. 1.1  Molecular  Diffusion 

The  first  quantitative  study  of  diffusion  was  made  by  Fiek  in  1855.  who  found 
an  analogy  between  molecular  diffusion  and  heat  conduction.  He  adapted  Fourier  s 
heat  equation  to  describe  the  diffusion  process  by  stating  that  mass  flux  of  a  diffusing 
substance  is  proportional  to  the  concentration  gradient, 

jo  =  -d0  VC,  (2.55) 


where  jo  is  the  mass  flux,  i.e.,  the  rate  of  mass  transfer  per  unit  area  of  cross  section, 
of  a  diffusive  non-reactive  tracer.  d0  is  the  coefficient  of  molecular  diffusion,  which  is 
a  scalar  in  aqueous  solution,  and  C  is  the  tracer  concentration,  usually  defined  as  the 
mass  of  tracer  per  unit  volume  of  solution.  In  the  absence  of  convection,  conservation 
of  mass  and  Equation  2.55  lead  to  the  well-known  Fickian  diffusion  equation  for 
transient  tracer  concentration. 

^  =  d0V2C,  •  (2.56) 

at 

where  t  is  time  and  V2  is  the  Laplacian  operator.  If  diffusion  occurs  within  a  fluid 
in  motion,  the  following  convection-diffusion  equation  governs  the  evolution  of  tracer 
concentration, 


dC  dC  „  , 

—  =  —  +  u  •  VC  =  d0V2C, 
dt  at 


(2.57) 


where  ^  -1-  u  •  V  is  the  Lagrangian  or  material  derivative  and  u  is  the  fluid 

velocity  vector. 

Albert  Einstein  in  1905  was  the  first  to  derive  an  expression  for  d0  of  a  colloidal 
suspension  in  an  aqueous  solution  (Einstein  1956).  His  derivation  was  based  on  the 
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assumptions  that  colloidal  molecules  are  approximately  spherical  and  are  large  in 
comparison  with  the  molecules  of  the  aqueous  solution.  The  successive  displacements 
of  the  colloidal  molecules  by  collision  with  liquid  molecules  were  independent  giving 
rise  to  a  random  walk  or  Brownian  motion.  Molecular  diffusion  (or  simply  diffusion) 
in  solids,  liquids,  and  gases  is  primarily  due  to  Brownian  motion.  It  is  also  caused 
by  osmotic  forces,  thermal  diffusion,  and  electro-osmosis.  Due  to  Brownian  motion, 
molecularly  dispersed  contaminants  tend  to  move  from  a  volume  element  with  a 
higher  concentration  or  specific  heat  toward  any  neighboring  element  with  a  lower 
concentration  or  specific  heat  content.  This  process  is  independent  of  convection. 

As  a  consequence  of  the  tortuous  diffusive  pathways  within  porous  media  and 
the  presence  of  fluid-solid  interfaces,  the  diffusion  coefficient  in  an  isotropic  porous 
medium  d  is  less  than  the  diffusion  coefficient  in  aqueous  solution  d0 , 

d  =  d0d*,  (2.58) 

where  d *  is  the  nondimensional  diffusivitv  of  the  medium.  From  a  review  of  data  on 
consolidated  granular  media  obtained  by  several  investigators,  Perkins  and  Johnston 
(1963)  suggested  that  the  value  of  d *  is  approximately  0.7  while  Fried  and  Combarnous 
(1971)  reported  values  of  0.4  to  0.8.  Bhattacharya  and  Gupta  (1990)  derived  a  formula 
for  d*  using  central  limit  theorem  applied  to  ergodic  Markov  processes,  which  provides 
a  theoretical  basis  for  the  values  of  d*  observed  in  experiments. 

Anisotropic  porous  media  have  different  diffusion  properties  in  different  directions. 
As  such,  the  direction  of  flow  of  diffusing  substance  at  any  point  in  such  media  may 
not  be  normal  to  the  surface  of  constant  concentration  through  the  point.  In  an 
anisotropic  medium,  the  vector  of  the  mass  flux  j  is  written  as, 

j  =  -dVC,  (2.59) 

where  the  diffusion  coefficient  d  is  a  symmetric  second-order  tensor  in  a  x-y-z  Carte¬ 
sian  coordinate  system, 
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dxx 

dxy 

dxz 

dXy 

dyy 

dy: 

d  x  2 

dyz  dZZ 

(2.60) 


Accordingly,  the  nondimensional  diffusivity  d*  is  also  a  symmetric  second-order  tensor 
in  the  coordinate  system. 


(2.61) 


Like  for  the  permeability  tensor,  rotation  of  the  coordinate  system  will  produce  a 
diagonal  d  tensor  when  the  coordinate  axes  are  aligned  with  the  principal  axes  of 
diffusion  for  the  medium.  Vectors  VC  and  j  are  parallel  in  such  cases  and. 


Ji  = 

h  = 
h  = 


,  dC 
_d,aT’ 

dC 

d 2  d-2 ’ 

,dC 

d 3  53’ 


(2.62) 

(2.63) 

(2.64) 


where  1,  2,  and  3  denote  the  three  principal  directions  of  the  diffusion  tensor,  ji  and 
di  (i  =  1,2,3)  are  the  mass  fluxes  and  diffusion  coefficients  along  the  three  principal 
directions,  respectively.  If  the  normal  of  a  surface  makes  angles  a,  /3,  and  7  with  the 
three  principal  axes  of  diffusion,  the  diffusion  coefficient  dn  in  the  direction  of  this 
normal  is  (Carslaw  and  Jaeger  1959), 


dn  =  d\  cos2  a  +  d2  cos2  (3  +  d3  cos2  7. 


(2.65) 


2. 3. 1.2  Mechanical  Dispersion 

Mechanical  dispersion  results  from  velocity  variations  within  the  porous  media. 
Three  mechanisms,  shown  schematically  in  Figure  2.9,  operate  at  the  microscopic 
level  to  give  rise  to  such  variations: 
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(1)  the  velocity  distribution  across  each  pore  channel; 

(2)  variations  in  velocity  as  a  result  of  the  distribution  of  pore  channel  size:  and 

(3)  tortuosity  of  individual  flow  paths. 

These  processes  combine  to  give  two  geometrical  aspects  of  dispersion:  (1 )  a  longitu¬ 
dinal  effect  due  to  the  differences  between  the  velocity  components  along  the  mean 
velocity  direction,  and  (2)  a  transverse  effect  due  to  the  differences  between  the  ve¬ 
locity  components  orthogonal  to  the  mean  velocity  direction.  Mechanical  dispersion 
is  attributed  to  convection.  When  convection  is  weak,  mechanical  dispersion  may  be 
negligible  relative  to  molecular  diffusion. 

In  most  dispersion  theories,  mechanical  dispersion  is  generally  considered  to  be 
mathematically  analogous  to  diffusion  because  spreading  of  the  contaminant  results 
from  velocity  variations  across  a  concentration  gradient.  Consequently,  the  mechani¬ 
cal  dispersion  component  of  tracer  flux  is  commonly  represented  mathematically  by 
a  ‘'Fickian  type"  equation  analogous  to  Equation  2.59, 

Jm  =  -DmVC,  (2.66) 

where  Jm  is  the  mass  flux  vector  due  to  mechanical  dispersion  and  Dm  is  the  coefficient 

of  mechanical  dispersion.  The  characterization  of  mechanical  dispersion  as  a  Fickian 
process  remains  a  subject  of  continuing  debate  which  will  be  discussed  in  detail  in 
section  2.3.7.  Adding  Equations  2.59  and  2.66  gives  the  total  dispersive  mass  flux 
vector  Jd  of  a  tracer  as, 

Jd  =  -D  VC,  (2.67) 

where  D  is  the  coefficient  of  hydrodynamic  dispersion  and, 

D  =  d  +  Dm.  (2.68) 

The  total  mass  flux  vector  J  due  to  both  convection  and  dispersion  is, 


J  =  uC  -  DVC. 


(2.69) 


Figure  2.9:  Components  of  mechanical  dispersion  at  the  microscopic  level:  (1)  the 
velocity  distribution  across  each  pore  channel;  (2)  variations  in  velocity  as  a  result  of 
the  distribution  of  pore  channel  size;  and  (3)  tortuosity  of  individual  flow  paths. 
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Conservation  of  mass  and  Equation  2.69  lead  to  the  widely  used  convection-dispersion 
equation  for  solute  transport  (Ogata  1970;  Bear  1972;  Freeze  and  Cherry  1979). 

d£  =  ?£  +  u.VC=DV2C.  (2.70) 

dt  dt 

The  convection-dispersion  equation  models  the  process  of  dispersion  as  a  diffusional 
process  and  predicts  that,  in  an  uniform  flow  field,  a  set  of  tracer  particles  will  be 
normally  distributed  about  a  center  moving  with  the  average  convective  velocity. 

From  the  viewpoint  of  classical  statistics,  if  the  travel  time  for  an  individual 
tracer  particle  becomes  much  larger  than  the  time  interval  during  which  its  successive 
velocities  are  positively  correlated,  its  total  displacement  may  be  considered  as  a  sum 
of  a  large  number  of  elementary  displacements  which  are  statistically  independent  of 
one  another.  Then,  the  probability  distribution  of  the  particle’s  total  displacement 
should  be  normal  according  to  the  central  limit  theorem  (Bear  1972).  In  view  of  the 
ergodic  principle  (in  an  ergodic  system,  the  ensemble  averages  and  time  averages  are 
equivalent),  this  distribution  also  represents  the  spatial  distribution  of  displacements 
of  a  cloud  of  initially  close  particles.  It  is  this  tendency  for  the  cloud  of  tracer  particles 
to  converge  to  a  normal  distribution  and  to  spread  with  a  variance  that  is  proportional 
to  time  that  makes  it  possible  to  model  mass  transport  as  a  diffusional  process. 

2.3.2  Hydrodynamic  Dispersion  Tensor  and  Its  Determination 

The  coefficient  of  hydrodynamic  dispersion  D  (Equation  2.67)  is  a  second-rank 
symmetric  tensor  (Bear  1972;  Koch  and  Brady  1985;  Plumb  and  Whitaker  1988).  The 
principal  axes  of  D  are  believed  to  be  oriented  parallel  and  transverse  to  the  mean 
direction  of  theTegional  flow.  This  indicates  that  mass  transport  can  be  defined  by 
two  characteristic  dispersion  components  that  are  defined  once  the  mean  direction  of 
regional  flow  is  known.  If  the  dispersion  coefficient  in  the  direction  of  the  flow  L  (i.e., 
the  longitudinal  dispersion  coefficient)  and  the  dispersion  coefficient  in  the  direction 
perpendicular  to  the  flow  T  (i.e.,  the  transverse  dispersion  coefficient)  are  denoted  as 
DL  and  Dt,  respectively,  the  convection-dispersion  equation  becomes, 
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?£  +  u  .  VC  =  Dl^j  +  DtV7tC, 

Ot  OLz 


(2.71) 


where  is  the  Laplacian  in  transverse  directions. 

Bear  (1961)  studied  the  relationship  between  D  and  u.  After  assuming  that  only 
part  of  each  velocity  component  is  significant,  which  is  either  parallel  or  normal  to 
the  mean  flow  direction,  he  showed. 


^  km  Q-ijki 


UjUj 

u0 


(t,  j,  k,  m.  —  x.  y.  z). 


(2.72) 


where  u$  is  the  magnitude  of  the  fluid  velocity,  ut  is  the  component  of  the  fluid 
velocity  in  i  direction,  and  is  called  the  media  s  geometrical  dispersivity  (a 

fourth-rank  symmetric  tensor).  In  general,  a^km  contains  81  components.  Seheidegger 
(1961)  studied  the  symmetry  properties  of  a^m  and  wrote  the  following  equation  for 
isotropic  media, 

aL  ~  aT  /r  r  ,  r  x  \  ^9 


d'ijkm  dT^ij^km  d” 


'  d"  8lm8jk  )  i 


where  ai  and  dx  are  the  longitudinal  and  transverse  dispersivities  of  the  media, 
respectively,  and  8  is  the  Dirac  delta  function. 


1  if  i  —  j\ 

0  if  i  ^  j- 

By  substituting  Equation  2.73  into  2.72,  the  following  equation  is  obtained, 

UiUj 


(2.74) 


Dij  —  ajUoSij  +  (ai  —  a?)- 


u0 


(2.75) 


which  provides  a  theoretical  basis  for  the  relation  between  D  and  u  for  isotropic 
media. 

In  practice,  the  longitudinal  and  transverse  mechanical  dispersion  coefficients  are 
generally  expressed  in  terms  of  the  seepage  velocity  vs  in  the  flow  direction  L, 

(2.76) 


DmL  =  ( aLvs)m\ 


DviT  —  i^T^s) 


m  2 


(2.77) 
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where  mi  and  m2  are  empirical  constants  between  1  and  2.  Laboratory  studies  have 
indicated  that  for  practical  purposes,  mi  and  m2  can  be  generally  taken  as  unity  for 
granular  materials  (Bedient  et  al.  1994).  For  isotropic  media,  from  Equations  2.68. 
2.58.  2.76,  and  2.77, 


Di  —  dod*  +  {aivs)m' ,  (2.iS) 

Dj  —  d(jd*  +  (ar^s)™2  .  (2.i  9) 

2 

As  most  chemical  species  have  do  on  the  order  of  10-9  to  10_1°  at  20oC.the 
contribution  of  diffusion  to  hydrodynamic  dispersion  is  typically  very  small  and  is 
neglected  in  most  models  of  groundwater  contamination.  One  exception  would  be  in 
geologic  formations  where  convective  transport  is  essentially  zero,  such  as  in  very  low 
permeability  clays  or  very  deep  aquifers  where  the  flow  rate  may  be  on  the  order  of 
a  few  centimeters  in  a  hundred  years. 

Dl  and  DT  are  also  defined  in  terms  of  the  spreading  of  an  initial  plume  (Fetter 


1993), 


(2.80) 

(2.S1) 


where  a\  and  a\  are  the  variance  of  the  longitudinal  and  transverse  spreading  of  the 
plume,  respectively.  In  the  asymptotic  limit  ( t  — >  oc),  Dl  in  Equation  2.80  and  Dr 
in  Equation  2.81  are  equivalent  to  those  defined  in  Equation  2.71  (Brenner  1980a). 

Many  experimental  studies  have  been  performed  to  evaluate  the  influence  on 
dispersion  of  the  characteristics  of  miscible  fluids,  of  the  flow  field,  and  of  media 
properties.  Conclusions  have  been  drawn  primarily  based  on  comparisons  between 
experimental  results  and  solutions  of  Equation  2.71.  The  studies  usually  lead  to  re¬ 
lationships  between  Dl,  Dt,  and  the  Peclet.  number  Pe.  Dimensional  analysis  shows 
that  D  is  a  function  of  Pe ,  which  defines  the  ratio  between  rate  of  transport  by 
convection  and  the  rate  of  transport  by  molecular  diffusion, 


(2.82) 
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Fried  and  Combarnous  (1971),  Bear  (1972),  and  Fried  (1975)  have  summarized  many 
of  the  experimental  results.  The  available  data  show  that  dispersion  phenomena 
exists  in  several  distinguishable  domains  according  to  the  value  of  Pe  (Figure  2.10). 
When  Pe  is  very  small,  molecular  diffusion  predominates;  in  the  next  domain,  the 


Figure  2.10:  A  generic  relation  between  Di  and  Pe. 

effect  of  molecular  diffusion  and  mechanical  dispersion  are  of  the  same  order;  in  the 
third  domain,  while  the  longitudinal  molecular  diffusion  is  negligible,  the  transversal 
molecular  diffusion  tends  to  reduce  the  longitudinal  spreading;  in  the  fourth  domain, 
mechanical  dispersion  predominates;  and  when  Pe  is  very  large,  while  mechanical 
dispersion  also  dominates,  the  effects  of  inertia  and  turbulence  can  no  longer  be 
neglected.  A  regressional  analysis  usually  leads  to  the  following  relation  between  Di 
and  Pe , 

=  aPe /3,  (2-83) 

do 

where  a  and  /?  are  empirical  constants  that  change  with  Pe. 

In  a  comprehensive  study  of  dispersion  phenomenon,  Klot.z  and  Moser  (1974) 
performed  2500  laboratory  column  tests  using  radioactive,  salt,  and  dye  tracers  to 
determine  the  dependence  of  Di  on  the  characteristics  of  the  permeant  and  porous 
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media.  Granular  materials  were  tested  with  differing  void  ratio,  uniformity  coefficient, 
and  grain  properties  (size,  shape,  and  angularity).  Their  study  generally  confirms  an 
approximate  linear  relation  between  Di  and  vs  as  suggested  by  Equation  2. 78  with 
nii  —  1-  In  addition,  the  effective  grain  size  and  uniformity  coefficient  substantially 
affect  the  value  of  Di.  The  influence  of  the  grain  shape  and  the  grain  roughness  as 
well  as  the  water  temperature  are  of  minor  importance.  In  the  study  on  transverse 
dispersion  in  column  tests  by  Klotz  et  al.  (1980).  they  evaluated  the  width  of  the 
tracer  cloud  by  measuring  the  concentration-time  distributions  perpendicular  to  the 
flow  direction.  In  general,  the  transverse  dispersion  was  difficult  to  determine.  They 
found  no  dependence  of  the  width  on  the  flow  velocity  within  the  low  measuring 
accuracy.  On  the  other  hand,  they  found  the  width  increased  with  increasing  grain 
size  of  the  porous  media.  The  tests  confirmed  that  Dt  is  6-20  times  smaller  than  DL 
and  anisotropic  mixing  occurs  even  in  isotropic  porous  media. 

Dispersion  theory  is  concerned  with  the  asymptotic  time-dependent,  spatial  dis¬ 
tribution  of  non-reactive  (passive)  tracer  particles  introduced  into  a  fluid  flowing 
through  the  interstices  of  a  porous  medium  under  the  influence  of  an  externally  ap¬ 
plied  pressure  gradient.  Research  on  dispersion  has  covered  numerous  aspects  and 
involved  many  approaches.  Various  models  have  been  developed  to  provide  a  math¬ 
ematical  description  of  this  process.  These  models  can  be  classified  into  four  groups: 
geometric,  random  network,  statistically  geometric,  and  statistical  models. 

2.3.3  Geometric  Model 

Taylor  dispersion  problem  is  probably  the  best  example  of  a  geometric  model. 
Taylor  (1953,  1954)  considered  the  dispersion  in  a  cylindrical  straight  tube  (Figure 
2.11).  Laminar,  steady,  one-dimensional  flow  of  a  fluid  (i.e.,  Poiseuille  flow)  takes 
place  in  the  tube.  The  velocity  distribution  is  parabolic  across  the  tube  with  an 
average  velocity  of  vs.  At  t  =  0,  a  tracer  is  introduced  at  the  origin  x  =  0,  simulating 
a  pulse  of  constant  concentration  C0.  The  tracer  particles  are  then  carried  away 
by  the  flow.  Because  of  the  parabolic  velocity  profile,  the  tracer  particles  dose  to 
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Figure  2.11:  Taylor  dispersion  problem. 


the  wall  move  slower  than  those  in  the  middle  of  the  tube.  The  spreading  of  the 
concentration  variation  front  is  purely  due  to  convection  initially.  However,  as  the 
front  moves  downward,  the  distribution  of  tracer  becomes  more  uniform  because 
of  transverse  molecular  diffusion  across  the  flow.  The  characteristic  time  for  this 
transverse  diffusion  is. 


(2.84) 


where  R  is  the  radius  of  the  tube.  Using  the  general  diffusion  equation  in  cylindrical 


coordinates,  Taylor  showed  that  for  t  >  tc,  the  average  cross  section  concentration 


C(x.t)  satisfies  the  following  equation, 


dC  dC  R2v2d2C 

dt  +  Vs  dx  0  +  48o?0  dx 2 


(2.85) 


Equation  2.85  indicates  that  longitudinal  dispersion  in  a  cylindrical  straight  tube 
obevs  the  same  law  as  molecular  diffusion  with  an  effective  dispersion  coefficient 

d2  2 

Drayior  —  do+  48^-  Tf*e  Taylor  dispersion  mechanism  is  one  in  which  the  transverse 
variation  of  longitudinal  velocity  and  transverse  mixing  interact  to  produce  an  overall 
longitudinal  mixing  process  that  is  Fickian.  Aris  (1956)  generalized  Taylor’s  approach 
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to  irregularly  shaped  tubes,  where  local  velocity  distributions  may  not  be  parabolic 
and  where  do  varies  with  concentration.  Using  a  moment-analysis  technique,  he 
obtained  the  coefficient  of  dispersion  D^ris  as. 


,  R2v; 

D.lris  —  do  +  a  — - — 


(2.86) 


where  R  is  a  characteristic  dimension  of  the  cross  section  and  a  is  a  dimensionless 
number  depending  on  the  cross  section. 

In  the  Tavlor-Aris  theory  of  dispersion,  the  dispersion  coefficient  is  proportional 
to  fj.  This  theory  cannot  be  applied  to  porous  media  in  general  because  miscible 
flow  in  porous  media  is  very  different  from  the  flow  in  a  capillary.  In  porous  media, 
most  of  dispersion  comes  from  the  meandering  of  streamtubes  through  the  complex 
void  structure,  and  not  from  the  presence  of  a  velocity  profile  within  each  pore. 

Brenner  (1980b,  1981,  1982)  coined  the  name  Taylor  dispersion  to  honor  Taylor’s 
work.  The  geometric  model  postulates  a  geometry  which  hopefully  bears  some  resem¬ 
blance  to  real  porous  media,  vet  is  sufficiently  simple  to  allow  analytical  treatment. 


2.3.4  Random  Network  Model 


Some  authors  have  developed  porous  media  dispersion  theory  by  representing  the 
media  as  random  networks  of  capillaries  (De  Josselin  de  Jong  1958;  Saffman  1959; 
Saffman  1960).  De  Josselin  de  Jong  (1958)  was  the  first  to  suggest  that  transverse 
dispersion  is  smaller  than  longitudinal  dispersion.  However,  his  assumptions  are  very 
restrictive:  he  neglected  molecular  diffusion  and  assumed  that  fluid  velocity  is  con¬ 
stant  within  each  capillary.  Saffman’s  model  (Saffman  1959;  Saffman  1960)  was  more 
general  and  represented  a  major  step  in  the  description  of  hydrodynamic  disper¬ 
sion.  The  model  consists  of  a  network  of  randomly  oriented  and  distributed  straight 
pores  in  each  of  which  the  flow  is  uniform  (i.e.,  constant  fluid  velocity).  The  pores 
are  connected  with  one  another  at  the  ends,  and  several  pores  may  start  or  end  at 
each  of  these  junctions.  The  path  of  a  fluid  particle  may  be  regarded  as  a  random 
walk  in  which  the  length,  direction,  and  duration  of  each  step  are  random  variables. 
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Using  this  approach,  Saffman  calculated  a  probability  distribution  function  for  the 
displacement  of  a  particle  in  time  and  deduced  values  for  the  dispersion  coefficients. 
By  assuming  the  pores  had  equal  circular  cross  sections  of  radius  R.  the  poie^  were 
all  of  equal  length  L,  and  an  average  velocity  vs  in  the  pores,  he  obtained  the  coeffi¬ 
cient  of  longitudinal  dispersion  Dsaffman,L  and  the  coefficient  of  transverse  dispersion 


D Saffman. T  as. 

Dsaffm  an.L 


do  3  R2vs2  ,  Lzv 

¥ +  80  To  -r  4 


LW  /\,  2 
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,  M  c-oth  M  -  1 
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d.v. 


(2  .Si 


and. 


d0  1  R2vs2  9 L2vs2  fl 

Ds.^m,T  -  J  +  + 
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J  0 


,,  M  coth  AI  —  1  , 
C>  BSP - - 


(2.SS) 


q  T  i  3 R2v2x2 

respectively,  where  M  -  and  B  =  d0  4 - 16Jq  ■  Saffman's  model  was  the 

first  that  showed  the  existence  of  several  domains  of  dispersion  with  the  change  of  the 
Peclet  number  Pe  and  the  difference  between  transverse  and  longitudinal  dispersion 


effects. 

The  models  of  De  Josselin  de  Jong  and  Saffman  were  extended  by  Haring  and 
Greenkorn  (1970)  to  the  case  of  non-uniform  media  by  the  use  of  beta  distributions 
for  the  distributions  of  pore  radius  and  pore  length.  Capillary  pressure,  permeability, 
and  longitudinal  and  transverse  dispersion  coefficients  were  calculated  in  terms  of  the 
parameters  of  the  beta  distribution. 

2.3.5  Statistically  Geometric  Model 


The  local  volume  averaging  method  produces  a  statistically  geometric  model.  This 
approach  is  based  on  the  concept  of  representative  elementary  volume  (REV)  (Bear 
1972)  and  the  Slattery- Whitaker  averaging  theorem  (Slattery  1969;  Whitaker  1969). 

An  REV  in  porous  media  is  the  smallest  volume  that  yields  the  local  average 
properties  (e.g.,  porosity)  such  that  addition  of  surrounding  media  to  this  volume 
does  not  change  these  average  values  (Figure  2.12).  On  the  other  hand,  it  must  be 
sufficiently  larger  than  the  size  of  a  single  pore  that  it  includes  a  sufficient  number  of 
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pores  to  permit  the  computation  of  meaningful  statistical  averages  required  for  the 
continuum  concept. 


Figure  2.12:  Porosity  n  versus  its  averaging  volume  V  in  a  porous  medium. 


The  Slattery- Whitaker  averaging  theorem  relates  a  transport  equation  at  the  mi¬ 
croscopic  level  to  its  macroscopic  counterpart.  The  theorem  states  the  relationship 
between  gradients  of  averages  and  averages  of  gradients  for  functions  defined  in  both 
solid  and  void  phases  which  exhibit  a  discontinuity  at  the  phase  interface  (Figure 
2.13), 


where. 


(Vi/>)  =  ^  f  uhdA, 

"  so/ id  —  flu  id 

W  =  ~  f  4>dV  =  ±  [  yd\'\ 

’  JV  V  JVfluid 


(2.89) 


(2.90) 


(  )  refers  to  an  average  quantity,  ip  is  a  scalar,  vector,  or  second-order  tensor  asso¬ 

ciated  with  the  fluid,  V  denotes  the  averaging  volume,  Vsoud  and  Vjiua  are  the  solid 
and  fluid  volumes  within  V,  respectively,  A  is  the  averaging  area  that  encloses  the 
volume  V,  Asoud-  fluid  is  the  solid-fluid  interfacial  area  within  .4,  and  n  is  a  unit  vector 
normal  to  Asoiid-.jluid.  Equations  2.89  and  2.90  make  it  possible  to  locally  average  the 
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Figure. 2. 13:  Representative  elementary  volume  (REX'). 

Navier-Stokes  equations  and  the  convection-diffusion  equation  on  an  REV.  A  detailed 
derivation  of  this  approach  can  be  found  in  Whitaker  (1969.  1973.  1986).  Bear  (19  <  2), 
and  Gray  (1975). 

By  locally  averaging  quantities  in  the  REX'*,  a  continuum  model  is  established  on 
the  scale  of  the  REV.  This  model  is  an  assemblage  of  randomly  interconnected  chan¬ 
nels  of  varying  length,  cross  section,  and  orientation.  However,  due  to  the  complexity 
of  the  flow  paths  and  the  interpore  and  pore-to-pore  fluid  dynamic  interactions,  a 
large  number  of  empirical  factors  must  be  introduced  in  the  derivations.  Statistically 
geometric  models  ignore  the  microscopic  nature  of  the  dispersion  process  dispersion 
occurs  not  in  one  continuous  medium,  but  in  a  medium  which  exhibits  abrupt  changes 
in  fundamental  properties  of  its  different  constitutional  phases.  This  approach  is  gen¬ 
erally  regarded  as  robust,  but  the  principal  weakness  is  that  rigorous  proof  of  the 
existence  of  satisfactory  averaging  volumes  (obeying  certain  ergodic  and  invariance 
requirements)  exists  only  for  idealized  porous  media.  In  fact,  porous  media  exhibiting 
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heterogeneity  on  a  continuum  of  scales  have  no  REV  (Cushman  1990:  Cushman  et  al. 
1994). 

2.3.6  Statistical  Model 

Probability  theory  and  statistical  mechanics  have  been  long  used  to  study  mass 
transport  in  porous  media.  Bear  (1972)  argued  that  the  path  of  a  tracer  particle  can 
be  visualized  as  the  vector  combination  of  two  motions:  one  along  the  pathline  of  a 
liquid  particle,  and  another  between  pathlines  as  a  molecular  diffusion  process.  The 
nature  of  both  motions,  the  first  determined  by  the  intricate  internal  geometry  of  the 
medium,  and  the  second  by  the  random  character  of  molecular  diffusion,  prevents 
any  deterministic  prediction  of  the  paths  of  tracer  particles.  In  a  similar  context. 
Greenhorn  and  Kessler  (1969)  stated  that  a  deterministic  description  of  dispersion 
is  useless  and  hopeless  in  practice  for  at  least  two  reasons:  (1)  determination  of 
the  precise  fluid-solid  boundary  of  a  random  porous  medium  is.  and  will  remain, 
impossible;  and  (2)  the  tortuous  boundary  (even  if  known)  will  render  the  problem 
mathematically  intractable. 

The  basic  postulate  of  the  statistical  approach  is  that,  although  it  is  impossible 
to  predict  the  exact  path  of  an  individual  tracer  particle,  one  may  employ  probability 
theory  to  predict  the  spatial  distribution  at  any  time  of  a  cloud  of  tracer  particles 
that  are  initially  at  a  close  proximity,  and  that  move  under  same  average  condition 
(Bear  1972).  If  ergodic  hypothesis  applies,  the  problem  of  the  spreading  of  a  cloud 
of  labeled  particles  is  reduced  in  the  statistical  approach  to  the  problem  of  random 
motion  of  a  single  tracer  particle  through  an  ensemble  of  random  porous  media.  The 
simplest  statistical  model  of  dispersion  is  the  one-dimensional  random  walk  model. 
In  this  model,  a  particle  moves  along  a  straight  line  as  a  series  of  steps  of  equal 
length,  and  each  step  is  taken  either  in  the  forward  or  backward  direction  with  equal 
probability  of  It  can  be  shown  that  if  the  particle  undergoes  n  displacements  per 
unit  time,  the  probability  that  it  will  be  between  x  and  x  4-  Ax  at  time  t  is, 
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p<I’t)Al=(5vbyexp('J^))Al’  <2'91) 

where  D  =  and  L  is  the  length  of  each  step.  Using  the  law  of  large  numbers. 
Scheidegger  (1954)  showed  that,  for  miscible  tracers,  concentration  at  a  point  in  the 
displacing  fluid  is  equal  to  the  probability  of  finding  a  particle  of  displacing  fluid  at 
that  point.  He  also  extended  the  random  walk  theory  to  three  dimensions.  However, 
his  analysis  led  to  a  scalar  D  resulting  from  his  negligence  of  molecular  diffusion  and 
his  assumption  of  isotropic  spreading. 

Unlike  models  which  require  non-physical  restrictions  on  the  allowable  dynamics 
or  media  heterogeneity,  or  both,  Cushman  (1997)  developed  a  general  dispersion  the¬ 
ory  using  classical  statistical  mechanics  which  is  fully  rigorous  and  obtained  without 
any  approximation  within  the  conceptual  framework  of  Hamiltonian  dynamics.  The¬ 
ories  in  statistical  mechanics  provide  the  tools  to  go  from  the  atomistic  or  moleculai 
structure  of  matter  to  a  continuum  scale.  Such  theories  allow  detailed  molecular  data 
to  be  mapped  into  field  and  constitutive  properties  or  material  parameters  (Cushman 
1997).  Cushman’s  nonequilibrium  statistical  mechanical  theory  of  transport  involves 
both  diffusive  and  convective  mixing  (dispersion)  at  all  scales.  The  results  are  based 
on  a  generalization  of  classical  approaches  used  in  molecular  hydrodynamics  and  on 
time-correlation  functions  defined  in  terms  of  nonequilibrium  expectations.  The  re¬ 
sulting  constitutive  laws  are  nonlocal  and  constitutive  parameters  are  wave  vectoi 
and  frequency  dependent. 

2.3.7  Probleip  of  Scale  and  Fickian  Approximation 

The  convection-dispersion  transport  equation  (Equation  2.70)  forms  the  theoret¬ 
ical  basis  in  almost  every  analysis  of  mass  transport  in  groundwater  systems.  This 
approach  has  traditionally  considered  dispersivity  as  a  characteristic  single-valued 
property  of  the  entire  medium  (Bear  1972).  However,  several  studies  suggest  that 
dispersivity  is  not  a  constant  but  rather  depends  on  the  mean  travel  distance  and 
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scale  of  the  problem  (Fried  1975;  Lallemand-Barres  and  Peaudecerf  1978;  Pickens  and 
Grisak  1981:  Gelhar  et  al.  1992).  Field  experiments  have  revealed  that  dispersivity 
increases  with  distance  between  the  source  and  the  point  of  observation  as  the  tracer 
encounters  progressively  larger  heterogeneities  within  the  aquifer.  Lallemand-Barres 
and  Peaudecerf  (1978)  were  the  first  to  publish  dispersivity  values  as  a  function  of 
the  distance  traveled.  A  plot  of  cii  versus  distance  showed  the  following  relationship. 

“■  -  w  {-02) 

where  ,r  is  the  distance  traveled  by  the  contaminant.  Equation  2.92  is  known  as  the 
"one-tenth  rule"  and  has  been  recommended  by  the  EPA  as  a  means  to  estimate 
dispersivity  in  preliminary  studies  where  no  data  exists  (U.S.EPA  1985).  In  a  com¬ 
prehensive  review  of  data  from  59  field  studies,  Gelhar  et  al.  (1992)  found  that  there 
was  a  systematic  increase  of  aL  with  observation  scale.  However,  they  stated  that  it  is 
not  appropriate  to  represent  the  relationship  between  ai  and  x  with  a  single  universal 
line  because  of  the  differing  degrees  of  aquifer  heterogeneity  at  different  sites.  Most 
published  values  of  aquifer  dispersivities  were  obtained  bv  calibrating  a  mass  trans¬ 
port  model  to  field  data.  This  has  usually  resulted  in  relatively  large  values,  typically 
on  the  order  of  tens  of  meters  to  over  100  meters.  The  magnitude  of  the  dispersivity 
is  a  measure  of  uncertainty  regarding  heterogeneous  flow  phenomena  in  an  aquifer. 
The  larger  the  degree  of  heterogeneity,  the  larger  will  be  the  dispersivity  if  the  system 
is  treated  as  homogeneous  with  respect  to  velocity.  If  a  non-homogeneous  velocity 
distribution  due  to  large-scale  geologic  heterogeneities  is  determined,  dispersivities 
would  depend  on  small  pore-scale  variations  of  velocity,  resulting  in  values  on  the 
order  of  a  few  centimeters  or  less.  Dispersivities  in  this  range  are  typical  of  column 
experiments  in  which  there  are  no  large-scale  heterogeneities  (Gelhar  et  al.  1992). 

As  can  be  seen,  the  domain  of  research  in  the  field  of  dispersion  is  part  of  the 
problem  of  change  of  scale,  which  is  a  very  general  problem  in  physics  (Fried  1975). 
Recent  theoretical  developments  have  shown  that  multiple  space-time  scales  appear 
naturally  in  the  description  of  the  transport  of  miscible,  non-reactive  solutes  in  sat¬ 
urated  porous  media  (Cushman  1984;  Cushman  1986;  Koch  and  Brady  1987;  Bhat- 
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tacharva  and  Gupta  1990;  Cushman  1997).  Dispersion  phenomena  manifest  different 
behavior  depending  on  different  scales.  Bhattacharya  and  Gupta  (1983)  showed  that 
mechanical  dispersion,  pore  diffusion,  and  molecular  diffusion  are  manifested  on  three 
separate  scales.  Scales  of  interest  often  are  labeled  as  microscopic,  pore,  macroscopic, 
full  aquifer,  and  regional.  Associated  with  each  scale  there  is  a  separate  transport 
equation  which  can  be  obtained  by  filtering  of  a  lower  scale  equation  (Cushman  1954; 
Cushman  1986). 

The  dispersion  problem  is  also  a  measurement  problem.  The  concept  of  scale  is 
fundamentally  connected  with  the  measurement  process  (Cushman  1986).  In  prac¬ 
tice.  tracer  concentrations  are  measured  by  some  means.  One  should  first  decide  on 
the  scale  of  interest;  and  then,  design  an  instrument  to  measure  on  that  scale.  The 
results  from  any  solution  of  the  transport  equations  are  valid  only  with  respect  to  the 
a  priori  scale  and  instrument  (Cushman  1986).  For  example,  laboratory  analysis  of  a 
tracer  sample  leads  to  a  concentration  value  averaged  over  a  scale  determined  by  the 
.size  of  the  sampling  device.  Baveye  and  Sposito  (1984)  attached  an  operational  signif¬ 
icance  to  the  averaging  associated  with  measurement  devices  and  suggested  defining 
a  macroscopic  variable  by  averaging  microscopic  values  according  to  a  weighting  func¬ 
tion  associated  with  the  sampling  device.  Moltyaner  (1989)  determined  the  weighting 
function  associated  with  the  so-called  through-the-wall  concentration  measurement 
device  by  considering  the  interaction  of  gamma  radiation  with  aquifer  materials  in  a 
laboratory  experiment.  The  function  was  found  to  follow  an  exponential  attenuation 
law. 

Some  approaches  have  been  proposed  to  model  scale-dependent  dispersion  (Pick¬ 
ens  and  Grisak  J.981;  Tompson  1988).  In  the  work  of  Pickens  and  Grisak  (1981), 
dispersivity  varies  temporally  as  a  function  of  mean  travel  distance  and  approaches 
a  maximum  or  asymptotic  value.  Tompson  (1988)  described  a  second-order  rela¬ 
tionship  for  local  dispersive  transport  which  can  be  cast  in  the  form  of  a  standard 
Fickian  relationship  with  apparent  time-dependent  dispersivity  functions  that  grow 
to  finite,  asymptotic  values.  To  model  scale-dependent  dispersion,  the  conventional 
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practice  is  to  involve  a  simple  scaling-up  of  the  porous  medium  dispersivity  value 
from  very  small  values  observed  in  column  experiments  to  the  much  larger  values 
calculated  from  field  trials.  However,  some  literature  suggest  that  dispersivity  may 
never  approach  a  finite  value  asymptotically  for  some  cases.  As  such,  the  Fiekian 
approximation  of  dispersion  and  the  traditional  convection-dispersion  equation  are 
not  valid  for  these  cases.  Matheron  and  Marsily  (1980)  showed  for  the  special  case 
of  a  stratified  porous  medium  with  flow  parallel  to  the  bedding,  the  transport  of  so¬ 
lute  cannot,  in  general,  be  represented  by  the  usual  convection-dispersion  equation, 
even  for  large  time.  Smith  and  Schwartz  (1980)  concluded  that  field-scale  dispersion 
(macrodispersion).  which  is  usually  caused  by  mixing  due  to  spatial  heterogeneities 
in  the  permeability  field,  cannot  be  modeled  by  a  large  value  diffusion  process.  Koch 
and  Brady  (1987)  showed  that  when  the  length  and  time  scales  on  which  a  transport 
process  occurs  are  not  much  larger  than  the  scales  of  variations  in  the  velocity  field 
experienced  by  a  tracer  particle,  a  description  of  the  transport  in  terms  of  a  local, 
average  Fiekian  process  is  not  applicable.  In  fact,  in  Cushman's  nonequilibrium  non¬ 
local  theory  of  transport  (Cushman  1997),  the  concept  of  classical  Fiekian  dispersive 
flux  is  a  special  case. 

2.4  Transport  Through  Spatially  Periodic  Porous  Media 

The  analysis  of  transport  processes  in  spatially  periodic  porous  media  is  a  rel¬ 
atively  simple  problem  when  numerical  or  analytical  calculations  are  confined  to  a 
unit  cell.  However,  if  one  were  to  analyze  an  irregular  unit  cell  of  arbitrary  shape, 
the  analysis  would  be  no  easier  than  that  of  other  models.  The  simplest  spatially 
periodic  porousjnedium  consists  of  a  two-dimensional  array  of  circular  cylinders  or 
three-dimensional  array  of  spheres.  Despite  its  simplicity,  no  rigorous  solutions  to  the 
transport  problem  are  available  for  irregular  arrays.  However,  some  are  available  for 
regular  arrays  of  cylinders  and  spheres. 
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2.4.1  Single-Phase  Flow 

The  flow  problem  is  the  steady  state  creeping  flow  of  an  incompressible  viscous 
Newtonian  fluid  through  a  rigid  spatially  periodic  poious  medium.  The  equations  to 
be  solved  are  the  Stokes  equations. 

0  =  -VP  +  /iV2u:  (2-93) 

V  •  u  =  0.  (2-94) 

These  equations  should  be  supplemented  by  the  no-slip  boundary  condition  at  the 
surface  Sp  of  the  solid  particles. 

u  =  0  on  Sp.  (2.95) 

In  addition,  another  fundamental  hypothesis  is  needed,  i.e.,  the  local  fluid  velocity  u 
is  spatially  periodic  with  the  same  periodicity  as  the  porous  medium. 

u(R)  =  u(R  +  Rn).  (2-96) 

This  hypothesis  is  based  on:  (1)  as  the  pressure  gradient  is  constant  across  the  unit 
cell,  it  can  be  expected  that  the  flux  is  also  constant;  and  (2)  locally,  u  must  satisfy 
the  no-slip  boundary  conditions  which  is  inherently  spatially  periodic. 

Hasimoto  (1959)  was  the  first  to  calculate  successfully  periodic  fundamental  so¬ 
lutions  of  the  Stokes  equations  for  the  flow  past  a  periodic  array  of  obstacles.  His 
approach  is  based  on  Fourier  series  expansion  and  he  applied  the  solution  to  dilute 
cubic  arrays  of  spheres  and  square  arrays  of  cylinders.  Happel  (1959)  employed  a  free- 
surface  model  to  study  flows  parallel  and  perpendicular  to  square  arrays  of  cylindeis. 
Snyder  and  Stewart  (1966)  calculated  the  velocity  and  pressure  profiles  for  Newtonian 
creeping  flow  through  dense  cubic  and  simple  cubic  beds  of  spheres  using  Galerkin  s 
method.  Their  solutions  converged  slowly  and  did  not  satisfy  continuity  exactly. 
Sorensen  and  Stewart  (1974)  improved  their  method  by  using  a  three-dimensional 
set  of  stream  functions  and  a  variational  principle.  Zick  and  Homsy  (1982)  started 
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from  Hasimoto's  fundamental  solution  and  obtained  a  system  of  Fredholm  integral 
equations  of  the  first  kind  for  the  unknown  stress  vector  at  the  surface  of  a  sphere, 
and  were  able  to  lower  the  dimension  of  the  problem  as  a  result.  Sangani  and  Acrivos 
(19S2b)  extended  the  analytical  technique  of  Hasimoto  and  were  able  to  deduce  a 
scheme  that  was  valid  for  the  whole  range  of  solid  concentrations  of  cubic  arrays  of 
spheres.  They  also  studied  the  cases  of  two-dimensional  square  and  hexagonal  arrays 
of  cylinders.  Drummond  and  Tahir  (1984)  adapted  the  method  of  singularities  to 
biharmonic  equations  to  study  Hows  parallel  and  perpendicular  to  arrays  of  cylinders. 

Due  to  computational  limitations  and/or  the  drawbacks  of  the  approaches  them¬ 
selves.  most  numerical  approaches  are  limited  to  study  two-dimensional  problems. 
Sangani  and  Acrivos  (1982a)  developed  a  numerical  technique  to  study  creeping  flow 
through  square  and  hexagonal  arrays  which  is  suitable  for  the  whole  range  of  poros¬ 
ity.  Larson  and  Higdon  (1986,  1987)  employed  the  boundary-integral  method  to  study 
microscopic  axial  and  transverse  flow  near  the  surface  of  media  composed  of  circular 
and  elliptical  inclusions.  Sangani  and  Yao  (1988a)  developed  a  numerical  method  by 
extending  Hasimoto's  approach  for  evaluating  macroscopic  transport  coefficients  by 
computing  the  many-particle  interactions  in  systems  with  an  arbitrary  size  and  spatial 
distribution  of  cylinders.  They  were  able  to  provide  creeping  flow  permeability  even 
for  some  random  spatially  periodic  structures.  Using  the  finite  element  method,  Mee- 
goda  et  al.  (1989)  studied  the  effect  of  specific  surface  area,  void  ratio,  particle  shape, 
material  heterogeneity  and  arrangement  of  particles  on  the  permeability  of  granular 
media.  They  proposed  an  equation,  which  is  similar  to  the  Kozeny-Carman  equa¬ 
tion.  to  predict  permeability.  Edwards  et  al.  (1990)  also  employed  the  finite  element 
technique  to  study  flow  fields  within  spatially  periodic  arrays  of  cylinders  arranged  in 
square  and  hexagonal  lattices  with  microscale  Reynolds  numbers  ranging  between  0 
and  200.  Bruschke  and  Advani  (1993)  developed  a  closed-form  solution  numerically 
for  flow  across  a  periodic  array  of  cylinders  by  matching  the  analytic  solution  using 
the  lubrication  approach  for  low  porosities  and  the  analytic  cell  model  solution  for 
high  porosities.  Ghaddar  (1995)  determined  the  statistical  transverse  permeability  of 
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two-dimensional  random  arrays  of  cylinders  by  employing  nip-elements  and  a  paral¬ 
lel  finite  element  procedure.  Koch  and  Ladd  (1997)  employed  the  lattice-Boltzmann 
method  to  study  the  effects  of  fluid  inertia  on  the  pressure  drop  required  to  drive 
fluid  flow  through  periodic  and  random  arrays  of  cylinders. 

Some  of  the  calculations  and  major  results  of  abovementioned  work  are  briefly 
outlined  below.  To  solve  Equations  2.93  to  2.96  for  an  incompressible  flow  through  a 
periodic  array  of  obstacles  with  their  centers  at  Rn,  a  useful  approach  is  to  continue 
the  interstitial  fields  analytically  into  the  space  occupied  by  the  obstacles,  replacing 
the  obstacles  with  singular  multipole  force  distribution  (Hasimoto  1959).  Equation 
2.93  is  replaced  by, 

0  =  -  VP  +  pV2u  -  Fn  £  <5(R  -  Rn),  (2-97) 

n 

where  Fn  is  the  force  acting  on  one  of  the  obstacles,  and  S  denotes  Dirac  delta  function 
defined  as.  •  -  • 

r  1  when  Rn  G  volume; 

/  <5(R  —  Rn)cf3R  =  <  (2-98) 

*'volume  I  0  when  Rn  ^  volume, 

and. 

5(R  -  Rn)  =  0  for  R  ^  Rn-  (2.99) 

Hasimoto  (1959)  then  expanded  u  and  VP  in  Fourier  series, 

u  =  ^u^e-2"^'R)  (2.100) 

R,n 

(2.101) 


(2.102) 


VP  =  £VP*ne- 

Rn 


•27rz(Rn  -R) 


where. 


Rn  =  njli  +  n2l2  +  n3I3, 
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which  are  vectors  in  the  reciprocal  lattice  determined  by, 

I2  x  I3 


Ii  - 

h  = 

i3  = 


T  ”  1 

v  unit  cell 

l3xii 

T  ’ 

*  unit  cell 

ii  x  I2 


(2.103) 

(2.104) 

(2.105) 


1  unit  cell 

By  introducing  the  series  into  Equations  2.97  and  2.94,  the  periodic  fundamental 
solution  of  the  Stokes  equations  for  the  flow  past  a  periodic  array  of  obstacles  is 
obtained.  The  analysis  also  leads  to  a  relationship  that  the  drag  forces  acting  on  the 
obstacles  within  the  unit  cell  are  balanced  by  the  mean  pressure  gradient  driving  the 
fluid,  i.e., 


n  =  -VPK 


unit  cell • 


(2.106) 


Following  the  same  procedure  as  Hasimoto  (1959),  for  flow  through  simple  cubic 
arrays  of  spheres  with  an  average  velocity  of  vs  in  the  direction  of  x\.  Sangani  and 
Acrivos  (1982b)  presented  the  formal  solution  in  a  xi-x2-x3  Cartesian  coordinate 
svstem  as. 
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(2.108) 

(2.109) 

(2.110) 


dXi  Vunit  cell 

where  R  is  the  radius  of  spheres,  and  Fd  is  the  dimensionless  drag  force  of  the  sphere 
in  the  direction  of  X\  defined  as, 

Fd 


Fd  = 


67 xnvsR' 


(2.111) 
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where  Fd  is  the  drag  force  acting  on  the  sphere  in  the  xx  direction.  G,  H .  and  L  are 


differential  operators,  such  that, 


A/  (  -A-nm 


H  =  £  £  « 


A/=0  m= 0 


52n  /0\*»  /_£_ 

9xjn  \d£/  \  dr] 


{M  =  n  +  2m). 


12.1 121 


£  =  X2  +  lx3) 


rj  -  x2  -  tx3, 


(2.113) 

(2.114) 


and  the  unknown  coefficients  Anm,  Bnm,  and  Cnm  are  to  be  determined  by  applying 
the  no-slip  boundary  condition  at  the  surface  of  the  spheres.  Sx  and  S2  in  the  above 


equations  are  given  by, 
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unit  cell 


For  two-dimensional  flow  perpendicular  to  an  array  of  cylinders,  if  the  mean  flow 
has  an  average  velocity  of  vs  in  the  xx  direction,  the  velocity  components  are  (Sangani 


and  Acrivos  1982b), 


(2.117) 


where. 
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The  analysis  leads  to  results  about  permeability  or  dimensionless  drag  force  as 
functions  of  porosity  of  the  media.  For  flow  through  a  spatially  periodic  simple  cubic 
array  of  spheres,  the  permeability  k  is  related  to  the  dimensionless  drag  force  Fd  by. 


•  _  f  unit  cell 

6nRFd 

Hasimoto  (1959)  found  that. 


(2.120) 


Fd-1  =  1  -  1.7601o5  +  6-  1.5593d2 +  0(o5).  (2.121) 


where  6  —  1  —  n  is  the  solid  fraction  of  the  media.  Sangani  and  Acrivos  (19S2b) 
refined  the  above  result  as, 

f-1  =  1  -  1.760105  +  0  -  1.559302  +  3.979905  -  3.O7340T  +  O(o^).  (2.122) 


For  two-dimensional  flow  perpendicular  to  the  axes  of  circular  cylinders,  the  di¬ 
mensionless  drag  force  is  usually  redefined  as  in  the  literature. 


Fd 

AiijiVs  ’ 


(2.123) 


which  is  the  dimensionless  drag  force  per  unit  length  of  the  cylinder  and  the  perme¬ 
ability  k  is  related  to  Fd  by. 


\r 

i  _  '  unit  cell 

4tt Td  ' 

Sangani  and  Acrivos  (1982a)  found  that, 

Fj1  =  In  0-0-5  -  0.738  +  0  -  O.88702  +  2.O3803  +  O(04), 
for  square  arrays  and, 

Frf_1  =  In  0~0-5  -  0.745  +  0  -  O.2502  +  O(04), 


(2.124) 


(2.125) 


(2.126) 


for  hexagonal  arrays.  Drummond  and  Tahir  (1984)  proposed  that, 

F~l  =  In 0— 1  -  1.476  +  20  -  O.502  -  O.O5104  -  O.O77508  +  O(012) 


(2.127) 
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for  square  arrays, 

F"1  =  In  <p-1  -  1.498  +  20-  0.5 02  -  O.OO25106  +  O{0n).  (2.128) 

for  staggered  arrays,  and. 

F~ i  -  In  1  -  1.354  +  26  -  0.5<?2  -  o!35S03  -  0.812o6  +  0(o12).  (2.129) 

for  hexagonal  arrays.  Attention  is  needed  that  all  the  above  functional  forms  between 
FT 1  and  6  are  usually  only  valid  for  small  6'  be.,  high  porosity  media. 

a  ■' 

2.4.2  Hydrodynamic  Dispersion 

To  describe  tracer  dispersion  in  spatially  periodic  porous  media,  in  addition  to 
solving  Equations  2.93  to  2.96,  the  following  convection-diffusion  equation  is  needed. 

—  +  u  •  VC  =  d0V2C.  (2.130) 

dt 

Impenetrability  of  the  solids  to  solute  transport  through  their  surfaces  requires,  .  • 

n  •  (uC  —  d0 VC)  =  0  on  Sp,  (2.131) 

where  n  is  the  unit  normal  to  the  solid  surfaces. 

Brenner  (1980a)  and  Brenner  and  Adler  (1982)  established  a  theory  for  dispersion 
in  flow7  through  spatially  periodic  porous  media.  The  theory  is  based  on  the  concept 
of  Browmian  motion  of  particles  and  makes  use  of  the  method-of-moments  developed 
by  Aris  (1956)  and  later  extended  by  Horn  (1971).  The  theory  leads  to  a  general 
formula  for  the  dispersion  tensor.  Brenner  showed  the  tensor  should  be  calculated 
bv  first  solving  a  convection-diffusion  problem  in  the  periodic  cell,  i.e.,  the  B  field 
problem.  Brenner’s  theory  is  rigorous  in  that  no  ad  hoc  assumptions  are  needed 
except  for  the  assumption  of  periodicity  of  the  media.  The  approach  is  now  called  the 
generalized  Taylor  dispersion  theory  (Brenner  1980b;  Brenner  1981;  Brenner  1982). 
By  local  volume  averaging  on  the  periodic  cell,  Carbonell  and  Whitaker  (1983)  were 
able  to  calculate  the  dispersion  tensor  through  an  f  field  problem.  Koch  et  al.  (1989) 
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studied  the  effect  of  grain  periodicity  on  dispersion  phenomena  and  the  aspects  of 
dispersion  mechanisms  which  are  artifacts  of  the  periodicity  constraint.  They  showed 
that  for  a  square  array  of  cylinders  or  a  simple  cubic  array  of  spheres,  and  in  the  limit 
Pe  — >  oc,  Di  depends  quadratically  on  Pe  and  Dj  approaches  a  constant  value.  Mei 
(1992)  used  a  homogenization  method  to  study  dispersion  in  porous  media.  Starting 
from  the  Navier-Stokes  equations  for  the  fluid  and  the  convection-diffusion  equation 
for  the  tracer,  he  arrived  at  Darcy’s  law  and  the  convection-dispersion  equation. 

By  solving  the  f  field  problem  of  C'arbonell  and  Whitaker  (19S3).  Eidsath  et  al. 
(19S3)  carried  out  a  numerical  simulation  of  dispersion  in  the  flow  through  a  square 
array  of  cylinders.  Using  a  finite  element  method,  they  first  solved  for  the  flow  field, 
and  then  solved  the  convection-diffusion  equation.  Edwards  et  al.  (1991)  studied  the 
Taylor  dispersion  of  a  passive  solute  within  a  fluid  flowing  through  a  two-dimensional 
spatially  periodic  porous  media  using  finite  element  based  on  Brenner's  B  field  prob¬ 
lem.  The  effects  of  microscale  Peclet  and  Reynolds  numbers  on  the  longitudinal  and 
transverse  dispersivities  were  also  investigated.  Salles  et  al.  (1993)  examined  theoret¬ 
ically  dispersion  in  spatially  periodic  porous  media  and  presented  numerical  results 
based  on  Brenner’s  B  method  and  random  walk. 

Brenner’s  theory  (Brenner  1980a)  and  the  method  of  homogenization  (Mauri  1991; 
Mei  1992:  Mei  et  al.  1996)  for  dispersion  in  spatially  periodic  porous  media  are  briefly 
outlined  below.  Brenner’s  theory  is  based  on  Brownian  motion  theory  for  particles. 
Suppose  that  at  time  t  =  0,  a  Brownian  solute  is  introduced  into  a  steady  state  flow  at 
some  arbitrary  position  R'.  The  instantaneous  position  R  =  R(f|R')  of  the  particle 
is  a  stochastic  variable.  The  probability  density  that  the  particle  is  located  at  R  at 
time  t,  knowing^that  it  was  released  at  R'  at  time  0,  is  denoted  by, 

p(R,  £|R').  (2.132) 

If  the  impenetrability  of  the  solid  phase  applies,  the  solute  must  lie  somewhere  within 
the  fluid  volume  for  t  >  0,  i.e., 

pd? R  =  1  for  t  >  0, 

roo 

fluid 


(2.133) 


DO 


where  Vjfuid  denotes  the  fluid  volumes  in  all  unit  cells.  In  addition,  it  will  be  assumed 
that.. 

p-M)  as  |R  -  R'|  — t  oc,  (2.134) 

at  a  sufficiently  rapid  rate  with  distance  to  ensure  the  convergence  of  the  integral  in 
Equation  2.133.  Conservation  and  continuity  of  probability  density  requires  that  p 
satisfy  the  following  equation, 

^  +  V-J  =  <S(R-R')  ■*(*).  (2.135) 

ot 

where  J  is  the  flux  of  the  probability. 

J  =  up  —  d0Vp,  (2.136) 

and  satisfies. 

n  •  J  =  0  on  Sp.  (2.137) 

The  probability  density  will  evolve  according  to  Equations  2.133  to  2.137-  The  mo¬ 

ments  of  the  probability  density  are  written  as, 

M m  =  [  (R-R')mp(R,t|RV3R-  (2.13S) 

Jvp 

The  first  three  moments  are  physically  the  most  interesting.  Moreover,  the  transients 
that  occur  just  after  the  introduction  of  the  solute  in  the  fluid  phase  are  not  of  interest, 
but  rather  what  occurs  after  a  long  enough  time  so  that  the  tracer  particle  can  sample 
all  the  interstitial  space  of  the  unit  cells.  This  condition  is  expressed  as, 

»  1.  (2.139) 

‘-'D 

The  spatial  integration  indicated  by  Equation  2.138  can  be  decomposed  into  two 
steps,  integration  over  a  unit  cell  and  then  over  all  unit  cells.  Brenner  performed  the 
calculation  and  found  that, 


Mo  =  1. 


(2.140) 
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This  zeroth-order  moment  represents  the  total  probability.  The  first-order  moment 
is. 

Mi  =  vst  +  B  +  Exp,  (2.141'1 

where  Exp  denotes  terms  that  go  exponentially  to  zero  as  t  — »  oc  and. 

B  =  — - -  [  B(r)c?3r.  (2.142) 

'  unit  cell  JVunllee„ 

in  which  B(r)  is  a  field  quantity  which  satisfies, 

d0V2B  —  V  ■  (uB)  =  vs,  (2.143) 

and  is  subjected  to  the  boundary  condition, 

n  •  VB  =  0  on  Sp.  (2.144) 

Additionally.  B  must  satisfy  the  pair  of  jump  conditions  across  opposite  cell  faces. 

(B)  =  -<r>, 

(VB)  =  0, 

where  the  jump  in  value  of  the  general  field  F  is  defined  as, 

(F)  =  F(Rn,  r  +  It)  —  F(Rn,  r),  (2.147) 

and  I,  are  the  three  basic  lattice  vectors  (Figure  2.1).  The  first-order  moment  rep¬ 
resents  the  average  position  of  the  solute  after  a  long  time.  Brenner  also  showed 
that, 

^  g  2v5v^  +  v,B  +  Bvs  +  2  d°  f  VB*  •  VBd3r  +  Exp,  (2.148) 

dt  *  unit  celt  JVunit  cell 

where  the  superscript  *  denotes  a  post-transposition  operator  (i.e.,  A*^k  =  Aikj  in 
Cartesian  tensor  notation).  The  physical  meaning  of  the  second-order  moment  lies  in 
the  calculation  for  the  mean  square  displacement  dyadic  (R  —  R)2, 


(2.145) 

(2.146) 


(R-  R)2  =  M2  -  MiM,, 


(2.149) 


CM 


where  R  is  the  mean  displacement  of  the  solute.  Brenner  then  followed  Einstein's 
generalized  Brownian  motion  relation  for  eontinua  and  defined  the  dispersion  tensor 
as. 

D  s  -  — (R  —  R)2  =  f  VB*  •  VBd3r.  (2.150) 

2  (it  \  unit  cell  J Vunit  cel, 

Furthermore.  Brenner  (1980a)  rigorously  proved  that  the  dispersion  tensor  defined  in 
Equation  2.150  is  the  same  quantity  as  in  Equation  2.70  of  the  conventional  continuum 
dispersion  model  as  t  — *  oo. 

Mei  (1992)  rederived  the  key  results  of  Brenner  using  the  theory  of  homogeniza¬ 
tion.  which  is  claimed  to  be  a  rigorous  mathematical  procedure  particularly  suited  for 
periodic  materials.  It  is  based  on  two  assumptions:  (1)  there  are  two  vastly  different 
length  scales  in  the  media,  and  (2)  the  media  is  periodic.  The  microscopic  length 
scale  £m  and  the  macroscopic  length  scale  £  of  a  spatially  periodic  porous  medium 
are  easily  identified  as  the  two  different  length  scales  for  homogenization  analysis. 

r 

Using  the  ratio  of  the  two  length  scales  e  =  —jP-  C  1  as  the  small  ordering  parameter, 
the  perturbation  expansions  are  introduced, 


Ui  =  uf'1  +  +  e2rq-2)  H - > 

(2.151) 

P  =  P{0)  +  eP(1)  +  e2P(2)  +  •  •  •  , 

(2.152) 

C  =  C{0)  +eCw  +e2C{2)  +  ■■■  , 

(2.153) 

where  and  pb)  are  functions  of  x,-  and  x\  =  ex*,  and  depends  on  xt,  x',  ti  =  et, 
and  t2  =  e2t.  Substituting  the  expansions  into  the  corresponding  governing  equations 
and  boundary  conditions,  a  set  of  perturbation  equations  are  obtained  at  successive 
orders.  The  leading-order  equation  is  homogeneous;  either  the  solution  itself  or  the 
coefficient  of  the  homogeneous  solution  are  indeterminate  and  independent  of  £m. 
The  next-order  solution  enables  the  calculation  of  the  leading-order  equation  and  the 
constitutive  coefficients,  such  as  the  permeability  and  dispersion  tensors. 
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2.4.3  Spatially  Periodic  versus  Real  Porous  Media 

Spatially  periodic  porous  media  have  been  extensively  studied  for  two  reasons, 
one  is  the  mathematical  tractability  and  the  other  is  the  feasibility  of  numerical 
implementation.  Although  assuming  spatial  periodicity  of  porous  media  does  not 
stand  on  solid  physical  evidence,  the  advantage  of  this  deterministic  approach  (i.e.. 
by  specifying  geometry  of  solid  inclusions  in  the  unit  cell)  lies  in  the  achievement  of 
a  complete  and  detailed  understanding  of  the  functional  relations  existing  between 
the  microscale  and  macroscale  processes.  To  study  phenomena  occurring  in  spatially 
periodic  porous  media  not  only  has  academic  importance,  it  also  provides  insights 
for  real  heterogeneous  disordered  porous  media.  Actually,  the  degree  of  geometric 
complexity  that  one  is  able  to  include  in  the  attempt  to  simulate  real  porous  media 
is  limited  only  by  the  capacities  of  computers.  With  the  continuing  development 
of  computer  technology,  the  ability  to  model  real  porous  media  through  the  use  of 
spatially  periodic  models  is  approaching. 

The  drawback  of  this  approach  is  that,  by  assuming  spatial  periodicity,  a  discrete 
separation  of  scales  has  been  required,  i.e.,  periodicity  of  the  ‘•microstructure"  with 
a  “scale-of-observation"  much  larger  than  the  period  (Cushman  1997).  Owing  to  the 
costs  and  poor  technology  associated  with  sampling  natural  porous  media,  a  periodic 
unit  cell  is  not  known  to  exist  a  priori,  and  if  it  does  exist,  the  scale-of-observation  is 
not  known  such  that  the  unit  cell  appears  infinitely  small.  Then,  local  theories  break 
down  and  nonlocal  with  infinite  support  on  the  scale-of-observation  has  to  be  assumed 
for  the  processes  (Cushman  1997).  Periodicity  may  also  produce  some  artifacts  in 
transport  phenomena.  Koch  et  al.  (1989)  studied  the  different  mechanisms  that 
contribute  to  dispersion.  They  stated  that  at  high  Peclet  number,  the  mechanical 
dispersion  that  is  induced  by  the  stochastic  fluid  velocity  field  in  disordered  media 
and  is  independent  of  the  molecular  diffusion  is  absent  in  periodic  media  where  the 
velocity  field  is  deterministic. 
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In  general,  the  potential  applicability  of  some  specific  results  gleaned  from  the 
spatially  periodic  analysis  should  not  be  overlooked  (Brenner  1980a).  However,  con¬ 
clusions  drawn  concerning  transport  in  real  disordered  porous  media  from  studies  of 
ordered  porous  media  should  be  viewed  with  caution  (Koch  et  al.  19891. 

2.5  Computational  Fluid  Dynamics  Techniques 

Once  a  porous  media  model  for  certain  phenomena  is  developed,  the  necessary 
computations  are  carried  out  analytically  or  numerically.  Lsing  advanced  computa¬ 
tional  fluid  dynamics  techniques,  numerical  simulations  of  porous  media  phenomena 
have  made  significant  progresses.  Conventional  numerical  methods  include  finite  ele¬ 
ment.  finite  difference,  and  boundary  integral  methods.  Their  application  in  porous 
media  transport  can  be  found  in  many  books  (Thomasset.  1981;  Wang  and  Ander¬ 
son  1982:  Liggett  and  Liu  1983:  Peyret  and  Taylor  1985).  Lattice-gas  and  lattice- 
Boltzmann  automata  are  two  new  methods  which  offer  efficient  simulations  in  complex 
geometries. 

2.5.1  Lattice-gas  Method 

Broadwell  (1964)  was  probably  the  first  to  develop  an  automata-type  model  for 
fluid  flow.  Velocity  was  a  discrete  variable  in  his  model,  but  space  and  time  were  both 
continuous.  Later,  Hardy  et  al.  (1976)  developed  a  completely  discrete  model  for  fluid 
flow  on  a  square  lattice  which  could  mimic  several  features  of  real  flow  problems. 

Lattice-gas  (LG)  models,  also  called  cellular-automata  models,  are  large  lattices 
in  which  each  site  can  be  in  one  of  several  discrete  states.  The  lattice  is  populated  by 
particles,  and  the  Boolean  variables  describing  the  state  of  the  system  indicate  the 
presence  or  absence  of  the  particles  at  the  intersections  of  the  lattice.  The  evolution  of 
the  system  is  governed  by  a  set  of  collision  rules  that  determine  how  the  particles  move 
in  the  lattice,  and  how  they  are  scattered  once  they  collide  with  each  other.  Those 
rules  are  chosen  to  conserve  mass  and  momentum  for  the  system.  Both  time  and  space 
are  discrete,  and  usually  site  connections  are  made  between  nearest  neighbors  only. 
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Frisch  et  al.  (1986,  1987)  and  Wolfram  (1986)  showed  that  in  order  for  the  discrete 
equations  to  reduce  to  the  Navier-Stokes  equations,  two-dimensional  simulations  have 
to  be  performed  using  a  triangular  lattice. 

Rothman  (1988)  used  lattice-gas  automata  to  model  flow  through  an  artificial 
two-dimensional  porous  medium.  He  showed  that  flow  rate  through  the  medium  was 
proportional  to  the  average  hydraulic  gradient  across  the  boundaries  in  accordance 
with  Darcy’s  law.  In  addition,  a  lattice  scale  must  be  chosen  such  that  the  linear  di¬ 
mension  of  any  void  space  is  at  least  twice  the  mean  free  path.  i.e..  the  average  distance 
traveled  before  particle  collision  occurs.  Brosa  and  Stauffer  (1989.  1991)  and  Kohring 
(1991)  also  used  two-dimensional  lattice-gas  models  to  simulate  the  flow  through  ran¬ 
domly  generated  porous  structures.  For  systems  with  up  to  22  million  sites.  Kohring 

A  A 

(1991)  found  permeability  varied  with  porosity  according  to  k  oc  for  n  >  0.7. 

He  also  noted  that  a  detailed  study  of  permeability  as  a  function  of  particle  size  and 
size  distribution  is  needed.  Rothman  (1990)  and  Gunstensen  and  Rothman  (1991) 
used  two-dimensional  lattice-gas  to  model  multiphase  flow  through  simple  pore  struc¬ 
tures.  Baudet  et  al.  (1989)  used  lattice-gas  to  simulate  two-dimensional  dispersion 
of  tracer  particles  for  flow  between  parallel  plates.  Simulation  results  agreed  very 
closely  with  analytical  solutions  under  laminar  flow  conditions. 

Difficulties  arise  when  lattice-gas  is  applied  to  three-dimensional  problems.  While 
a  triangular  lattice  contains  sufficient  symmetry  to  yield  a  macroseopically  isotropic 
two-dimensional  fluid,  no  regular  lattice  exists  that  yields  isotropy  in  three  dimen¬ 
sions.  d’Humieres  et  al.  (1986)  theorized  that  flow  could  be  modeled  in  a  four¬ 
dimensional  face-centered  hypercubic  (FCHC)  lattice.  Three-dimensional  results  are 
obtained  by  projecting  the  four-dimensional  data  down  to  a  three-dimensional  space. 
Frisch  et  al.  (1987)  discussed  the  approach  in  detail.  The  FCHC  model  has  been 
implemented  by  Rivet  et  al.  (1988)  for  the  simulation  of  three-dimensional  vortex 
behind  a  circular  plate  in  an  uniform  flow  field. 
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2.5.2  Lattice-Boltzmann  Method 

Lattice-gas  not  only  encounters  difficulties  in  three-dimensional  problem,  it  also 
suffers  from  statistical  noise.  The  lattice-Boltzmann  (LB)  method  has  been  devel¬ 
oped  to  overcome  some  of  the  drawbacks  of  the  lattice-gas  method.  In  the  spirit  of 
the  Boltzmann  equation  of  kinetic  theory,  the  lattice-Boltzmann  technique  models  a 
fluid  according  to  the  average  behavior  of  particles  on  a  lattice  rathei  than  as  the  dis¬ 
crete  particles  themselves.  The  Boolean  site  populations  of  the  lattice-gas  automata 
become  real  numbers  between  0  and  1  representing  their  average  value  through  time 
as  the  simulation  progresses.  Thus,  instead  of  following  each  particle  in  detail,  the 
model  follows  the  average  behavior  of  particles  at  a  given  site.  This  also  eliminates 
the  statistical  noise  associated  with  the  lattice-gas  approach.  Lattice-Boltzmann  is 
particularly  efficient  for  low  Re  simulations  (McNamara  and  Zanetti  19SS)  and  is 
unconditionally  stable  for  the  solution  of  the  Navier-Stokes  equations  (Frisch  1991). 

Using  the  lattice-Boltzmann  FCHC  method.  Succi  et  al.  (1989)  simulated  Darcy's 
law  for  fluid  flow  through  a  three-dimensional  random  medium  with  32s  =  32  <68 
sites.  Pores  were  no  smaller  in  cross  section  than  4x4  lattice  units  and  no-slip 
boundary  conditions  were  imposed  at  the  walls.  For  three  values  of  porosity,  flow 
rate  was  found  proportional  to  hydraulic  gradient  in  accordance  with  Darcv  s  law 
when  Re  <  5.  Furthermore,  Succi  et  al.  (1991)  demonstrated  that  the  Boltzmann 
technique  is  a  viable  tool  to  simulate  the  flow  regime  for  a  variety  of  hydraulics 
problems  ranging  from  laminar  to  turbulent  conditions.  Cancelliere  et  al.  (1990) 
simulated  flow  through  a  three-dimensional  FCHC  ‘‘penetrable  sphere  model  using 
lattice-Boltzmann  method  to  calculate  permeability  as  a  function  of  the  solid  fraction 
of  the  media.  The  simulated  data  fit  the  Kozeny-Carman  equation  well.  Gunstensen 
and  Rothman  (1992,  1993)  used  a  lattice-Boltzmann  model  for  the  simulation  of 
immiscible  flow  of  two  fluids  through  a  three-dimensional  porous  medium  consisting 
of  randomly  overlapping  spheres.  Ladd  (1994a,  1994b)  simulated  the  hydrodynamic 
interaction  of  particles  in  a  fluid  suspension  using  the  lattice-Boltzmann  method. 
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2.5.3  Smoothed  Particle  Hydrodynamics  (SPH) 

Smoothed  Particle  Hydrodynamics  (SPH)  can  also  be  used  to  investigate  pore- 
scale  transport  phenomenon  in  porous  media.  First  developed  for  astrophysical  appli¬ 
cations  (Lucy  1977;  Gingold  and  Monaghan  1977).  SPH  has  been  applied  successfully 
to  a  wide  range  of  problems.  SPH  is  a  fully  Lagrangian  computational  fluid  dynamics 
technique  in  which  the  numerical  solution  is  achieved  without  a  grid.  I  sing  this  ap¬ 
proach.  fluid  velocity,  pressure  and  tracer  distributions,  discharge  velocity,  and  fluid 
particle  pathlines  can  be  computed,  as  well  as  other  information  that  would  be  dif¬ 
ficult  or  impossible  to  observe  experimentally  or  with  other  numerical  approaches. 
SPH  has  a  number  of  advantages  over  competing  numerical  techniques.  Mobile  solid 
boundaries  are  difficult  to  incorporate  into  more  conventional  methods,  which  require 
either  continual  remeshing  of  the  domain  or  complicated  algorithmic  modifications. 
The  meshless  nature  of  SPH  simplifies  the  simulation  of  mobile  (Monaghan  1994)  or 
even  deformable  boundaries  (Libersky  and  Petschek  1990;  Benz  and  Asphaug  1995: 
Randles  and  Libersky  1996).  The  Lagrangian  nature  of  SPH  simplifies  the  inclusion  of 
extra  physical  effects  at  a  fluid-fluid  boundary.  For  example,  it  is  possible  to  simulate 
immiscible  fluids  with  SPH  (Morris  1999),  which  is  of  crucial  importance  to  modeling 
the  mobility  of  non-aqueous  phase  liquids  within  a  solid  matrix.  Most  methods  suffer 
from  an  increase  in  complexity  when  extended  to  three-dimensional  problems.  The 
SPH  algorithm  remains  essentially  unchanged  when  considering  one,  two,  or  three 
dimensions.  In  addition,  most  formulations  of  SPH  guarantee  local  conservation  of 
mass,  momentum,  and  energy.  This  is  not  typically  the  case  with  competing  methods, 
such  as  the  finite  element  method.  While  SPH  is  versatile,  errors  can  sometimes  be 
larger  than  those  obtained  using  grid-based  methods  tailored  for  specific  problems. 
Moreover,  SPH  can  be  computationally  expensive  for  certain  applications,  although 
at  comparable  resolutions,  the  computational  expense  of  SPH  is  comparable  with 
conventional  methods.  SPH  has  been  extended  to  model  solid  dynamics  problems  as 
well  (Libersky  and  Petschek  1990;  Benz  and  Asphaug  1995;  Randles  and  Libersky 
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1996)  and  it  has  only  recently  been  applied  to  low  Reynolds  number  incompressible 
flows  (Morris  et  al.  1997;  Zhu  et  al.  1997:  Zhu  et  al.  1999). 

2.6  Summary 

A  systematic  study  of  pore-scale  transport  through  spatially  periodic  porous  media 
using  a  true  Lagrangian  method  has  not  been  reported  in  the  literature.  Although  the 
flow  problem  has  been  studied  extensively  using  both  circular  and  elliptical  inclusions 
(Hasimoto  1959;  Happel  1959;  Snyder  and  Stewart  1966:  Sorensen  and  Stew  ait  19 1 4. 
Zick  and  Homsv  1982;  Sangani  and  Aerivos  1982b:  Sangani  and  Acrivos  1982a:  Larson 
and  Higdon  1986;  Larson  and  Higdon  1987:  Sangani  and  Yao  1988a:  Sangani  and 
Yao  1988b;  Drummond  and  Tahir  1984;  Meegoda  et  ah  1989;  Edwards  et  ah  1990: 
Bruschke  and  Advani  1993;  Ghaddar  1995;  Ranganathan  et  ah  1996;  Koch  and  Ladd 

1997) .  only  a  few  papers  have  discussed  the  problem  of  diffusion  and  hydrodynamic 
dispersion  for  circular  inclusions  (Eidsath  et  ah  1983:  Koch  et  ah  1989:  Edwards 
et  ah  1991;  Salles  et  ah  1993).  However,  most  of  the  solutions  were  achieved  by 
either  solving  Brenner  (1980a)  B  field  problem  or  Carbonell  and  Whitaker  (1983)  f 
field  problem  without  obtaining  the  evolution  and  distribution  of  concentration.  The 
following  chapters  present  a  study  of  pore-scale  transport  through  spatially  periodic 
porous  media  using  Smoothed  Particle  Hydrodynamics.  Starting  from  basic  fluid 
flow  and  diffusion  equations,  the  flow  and  concentration  fields  are  solved  to  yield 
values  of  permeability  and  coefficient  of  diffusion  and  hydrodynamic  dispersion.  Using 
this  Lagrangian  approach,  all  relevant  information  regarding  flow,  diffusion,  tracer 
convection  and  hydrodynamic  dispersion  can  be  obtained. 
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CHAPTER  3  DEVELOPMENT  AND  VERIFICATION  OF  A 
PORE-SCALE  FLOW  MODEL  USING  SMOOTHED  PARTICLE 

HYDRODYNAMICS 

Smoothed  Particle  Hydrodynamics  (SPH)  was  originally  developed  for  astrophys- 
ical  applications  to  model  compressible  fluids  at  high  Reynolds  number  (Lucy  1977: 
Gingold  and  Monaghan  1977).  SPH  is  well  suited  to  model  compressible  flows  be¬ 
cause.  in  SPH,  the  fluid  is  driven  by  local  density  fluctuations  at  the  particles.  Mon¬ 
aghan  (1994)  extended  SPH  to  inviscid  incompressible  flow  problems  involving  free 
surfaces  for  high  Reynolds  numbers  and  free-slip  boundary  conditions.  Efforts  have 
also  been  made  to  simulate  compressible  gases  with  Reynolds  numbers  as  low  as  5 
(Takeda  et  al.  1994).  The  extension  of  SPH  to  model  low  Reynolds  number  (Re  <  1) 
incompressible  flows  calls  for  modifications  of  the  standard  SPH  formalism  which  are 
discussed  by  Morris  et  al.  (1997),  Zhu  et  al.  (1997),  and  Zhu  et  al.  (1999).  In  this 
chapter,  an  overview  of  SPH  is  first  provided,  the  SPH  numerical  approach  to  model 
porous  media  flow  is  discussed  in  detail  and  verification  of  the  model  is  presented. 

3.1  Overview  of  SPH 

The  standard  approach  to  SPH  is  reviewed  by  Benz  (1990)  and  Monaghan  (1992). 
In  SPH.  a  compressible  fluid  is  represented  by  a  field  of  disordered  particles  (Figure 
3.1),  typically  of  fixed  mass,  which  follow  the  local  fluid  motion,  convect  contact 
discontinuities,  preserve  Galilean  invariance,  and  reduce  computational  diffusion  of 
various  fluid  properties  including  momentum.  The  equations  governing  the  evolution 
of  the  fluid  become  expressions  for  interparticle  forces  and  fluxes  when  written  in  SPH 
forms.  In  standard  SPH,  each  particle  carries  mass  m ,  density  p,  velocity  u,  and  other 
fluid  quantities  specific  to  a  given  problem.  The  particle  is  mathematically  treated 
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Figure  3.1:  Sphere  of  influence  for  SPH  particle  a. 
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as  an  interpolation  point  at  which  fluid  properties  are  computed  as  the  weighted  sum 
of  values  from  neighboring  particles.  To  illustrate,  consider  a  field  quantity  .4(r) 
expressed  by. 

.4(r)  =  J  A(t')S(t  —  T')dTr ,  (3.1) 

where  r  and  t'  are  position  vectors  and  5  is  the  Dirac  delta  function.  If  6  is  replaced 
with  an  interpolation  kernel  ir(r./?),  an  integral  interpolant  .4j(r)  of  the  function  is 
obtained. 

.4;(r)  =  J  .4(r')TT(r  -  r\  h)dr'.  (3.2) 

The  kernel  typically  takes  the  form, 

(3-3) 

where  a  is  the  number  of  dimensions  for  the  problem  and  h  is  the  smoothing  length,  a 
natural  length  scale  associated  with  the  method  of  SPH.  The  kernel  has  the  following 
properties, 

J  H'(r  —  r',  h)dv'  —  1,  (3.4) 

and. 


lim  IT(r  —  r',  h )  =  J(r  —  r').  (3-5) 

h  — >0 

For  numerical  work.  .4j(r)  is  approximated  by  a  summation  interpolant  -4s(r)  over 
the  particle  field, 


r)  =  —AbW(\r  -  r6|,  h), 
6 


(3.6) 


where  field  quantities  at  particle  b  are  denoted  by  subscript  The  quantity  ^  is 
the  inverse  of  the  number  density  at  particle  b  and  can  be  considered  as  the  volume 
of  fluid  associated  with  that  particle. 
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Using  the  above  concepts,  SPH  equations  governing  fluid  motion  can  be  obtained. 
For  example,  the  fluid  density  pa  at  particle  a.  may  be  evaluated  by  direct  particle 
mass  summation  as, 

Pa  =  Yl  ^3'7- 

b 

where. 

Wab  =  \V(rab.h).  (3.SI 

and. 

rab  =  rQ  -  rb.  (3.9) 

Other  expressions  for  derived  field  quantities  at  the  particles  are  obtained  by  sum¬ 
mations  involving  the  kernel  and/or  its  derivatives.  As  derivatives  can  be  obtained 
by  ordinary  differentiations  provided  that  II  (r,  h)  is  differentiable,  there  is  no  need 
for  a  grid.  For  example,  the  gradient  and  divergence  of  .4(r)  can  be  obtained  by, 

Vdj(r)  =  Y  —  Ai,VI'F(|r  -  r6|,  h),  (3.10) 

V  Pb 

and, 

V  •  A,( r)  =  ~Ab  •  Viydr  -  r"l’ h)'  (3-H) 

b  Pb 

respectively.  However,  in  practice,  it  is  usually  more  accurate  to  use, 

VA  =  -  (V(pA)  -  AVp)  =  -  Ymb(Ab  -  A«)VH'(|r  -  r6| ,h),  (3.12) 

and.  . 

V  •  A  =  -  (V  •  ( pA )  -  A  ■  Vp)  =  —  Y  mb(Ab  -  Aa)  •  VVF(|r  -  rb\,h). 

P  Pa  b 


(3.13) 
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3.2  Model  Establishment  for  Pore-Scale  Flow  Problem 

Groundwater  flow  is  generally  regarded  as  incompressible  since  the  bulk  fluid 
velocity  is  much  smaller  than  the  corresponding  speed  of  sound.  The  solution  to  an 
incompressible  flow  problem  is  achieved  by  solving  the  following  mass  and  momentum 
conservation  equations  throughout  the  flow  domain  subjected  to  proper  boundary 
conditions. 

V  • u  =  0,  (3.14) 

—  =  --Vp  +  uV2u  +  g,  (3.1  o) 

at  p 

where  v  is  the  fluid  kinematic  viscosity. 

While  SPH  is  well  suited  to  model  compressible  fluids,  the  required  modifications 
and  extensions  on  the  method  of  SPH  for  modeling  porous  media  incompressible  flow 
include  treatment  of  real  viscosity,  introduction  of  dynamic  pressure  concept,  choice 
of  the  equation  of  state  and  kernel  interpolation,  and  implementation  of  no-slip  flow 
boundary  conditions. 

3.2.1  Equation  of  State 

SPH  cannot  model  a  truly  incompressible  fluid  as  in  SPH.  fluid  pressure  is  an 
explicit  function  of  local  fluid  density  and  particle  motions  are  driven  by  local  density 
gradients.  In  SPH,  an  incompressible  flow  must  be  approximated  by  solving  for  the 
flow  of  a  nearly  incompressible,  or  quasi-incompressible,  fluid.  As  a  result,  to  close 
the  group  of  equations  used  for  describing  compressible  fluids,  an  equation  of  state  is 
required  in  the  form  of, 

P  =/(*>).  (3.16) 

Theoretically  the  actual  state  equation  for  the  fluid  modeled  could  be  used  but  this 
usually  results  in  a  prohibitively  small  time  step  for  numerical  stability  due  to  CFL 
condition  (Courant  et  al.  1928).  Therefore,  in  SPH,  the  incompressible  fluid  flow 


69 


has  to  be  modeled  as  a  quasi-incompressible  flow  through  a  quasi-incompressible 
equation  of  state.  But  the  choice  of  the  numerical  speed  of  sound  should  control  the 
density  variation  within  a  reasonable  range.  In  this  work,  the  chosen  sound  speed  is 
low  enough  to  be  practical,  and  yet  high  enough  to  limit  fluid  density  fluctuations 
to  about  1%.  A  similar  approach  has  been  used  in  grid-based  techniques  to  model 
steady  incompressible'  flows  (Chorin  1967:  Turkel  1987:  Tamamidis  et  al.  1996). 

Most  grid-based  techniques  model  the  flow  of  water  as  incompressible  since  the 
speed  of  sound  in  water  is  usually  very  large  compared  with  the  bulk  fluid  motions 
(i.e..  a  very  low  Mach  number).  Previous  applications  of  SPH  to  incompressible  flows 
of  water  (Monaghan  1994:  Monaghan  1995b)  have  modified  a  state  equation  suggested 
by  Batchelor  (1967)  which  describes  sound  waves  accurately. 


where  7  =  7  and  a  zero  subscript  denotes  reference  quantities.  The  choice  of  7  =  7 
in  Equation  3.17  causes  pressure  to  respond  strongly  to  variations  in  density.  Thus, 
perturbations  to  the  density  field  remain  small,  even  at  high  Reynolds  numbers. 
However,  as  the  density  fluctuations  increase,  small  errors  in  density  correspond  to 
increasing  larger  errors  in  pressure.  For  lower  Reynolds  number  porous  media  flows, 
more  accurate  pressure  estimates  are  obtained  using  SPH  if  7  is  taken  to  be  unity  as 
in  grid-based  approaches  (Chorin  1967;  Turkel  1987;  Tamamidis  et  al.  1996),  since 
errors  in  density  and  pressure  remain  proportional. 

In  previous  work  involving  incompressible  fluids,  the  subtraction  of  1  in  Equation 
3.17  was  introduced  to  remove  spurious  boundary  effects  at  a  free  surface.  It  is  well 
established  that -SPH  is  unstable  when  attractive  forces  act  between  particles  (Mor¬ 
ris  1994;  Swegle  et  al.  1995;  Balsara  1995;  Morris  1996b).  Consequently,  for  flow 
simulations  in  this  work,  this  subtraction  was  found  to  lead  to  numerical  instabilities 
in  regions  of  sustained  low  pressure.  Since  the  simulations  (and  many  other  applica¬ 
tions)  have  particles  filling  all  space  (section  3.3.1),  the  following  artificial  equation 
of  state  is  used. 


70 


P  =  c2P, 


(3.18) 


where  c.  is  the  speed  of  sound. 

The  sound  speed  must  be  chosen  carefully  to  ensure  both  an  efficient  and  accurate 
solution  of  a  given  problem.  The  value  of  c  must  be  large  enough  that  the  behavior 
of  the  corresponding  quasi-incompressible  fluid  is  sufficiently  close  to  that  of  the  real 
fluid,  yet  it  should  not  be  so  large  as  to  make  the  time  step  prohibitively  small. 
Monaghan  (1994)  argued  that,  for  density  to  vary  by  at  most  lc/c.  the  Mach  number 
of  the  flow  which  is  defined  as. 


(3.19) 


should  be  0.1  or  less,  where  U  is  the  typical  fluid  velocity  scale  of  the  problem.  In 
fact,  for  typical  smoothing  lengths  used  with  SPH,  the  kernel  interpolation  is  only 
accurate  to  within  approximately  1%.  The  principal  cause  of  this  variation  is  small 
fluctuations  in  density  which  inevitably  occur  as  particles  move  past  one  another. 
Thus,  local  pressure  gradients  obtained  using  a  high  sound  speed  are  potentially 
noisy.  Nevertheless,  the  velocities  obtained  are  accurate  if  smoothed  either  by  XSPH 
(Monaghan  1989)  or  viscosity.  It  was  found  that  the  computed  pressure  field  is  in 
close  agreement  with  other  techniques  when  c  is  chosen  such  that  the  density  varies 
about  1%. 

Considering  the  balance  of  forces  in  Equation  3.15.  c2  should  be  comparable  with 
the  largest  of, 


U2  uU  FC 

T’  ex' 


(3.20) 


where  £  is  the  typical  length  scale  of  the  problem,  F  is  the  magnitude  of  the  driving 
body  force  of  the  flow,  and  A  is  defined  as, 

A  =  — , 

Po 


(3.21) 


where  A p  is  the  maximum  density  difference  in  the  flow  and  po  is  the  initial  fluid 
density.  A  value  of  A  =  1 7c  was  chosen  for  this  study  because  SPH  kernel  inter¬ 
polation  is  only  accurate  to  within  approximately  1%.  The  first  term  in  Equation 
3.20  corresponds  to  that  derived  by  Monaghan  (1994).  The  second  and  third  terms 
ensure  that  pressure  forces  are  comparable  with  viscous  and  body  forces,  respectively. 
Equation  3.20  provides  a  first  estimate  of  the  sound  speed  to  use  for  a  given  problem. 
To  obtain  a  better  estimate  of  c~.  an  initial  simulation  is  run  at  low  resolution  to  find 
the  distance  A L  in  the  direction  of  the  applied  body  force  between  the  locations  of 
minimum  and  maximum  fluid  density.  The  third  term  of  Equation  3.20  then  becomes. 

PoFAZ/  ^  22) 

Ap 

3.2.2  Evolution  of  Density 

If  Equation  3.7  is  used  to  evolve  density  when  modeling  incompressible  free  surface 
flows,  particle  density  is  smoothed  out  at  the  edge  of  the  fluid  and  spurious  pressure 
gradients  are  induced  at  the  surface.  To  avoid  this  problem,  Monaghan  (1994)  initially 
set  the  density  to  a  reference  value,  applied  a  more  accurate  estimate  of  the  divergence 
of  the  velocity  field,  and  evolved  particle  densities  according  to  the  following  SPH 
equation  for  continuity, 

^  =  £m6utt6-VaW'a6>  (3.23) 

b 

where  VQ  denotes  the  gradient  with  respect  to  the  coordinates  of  particle  a  and, 

ua6  =  ua  -  u6.  (3-24) 

Equation  3.23  is"derived  from  the  mass  conservation  equation  for  a  compressible  fluid, 

^  =  -pV  •  u.  (3.25) 

at 

The  simulations  of  interest  here  do  not  involve  free  surfaces  and  Equation  3.7 
may  be  used  to  evolve  particle  densities.  One  disadvantage  of  this  approach  is  that 
density  must  be  evaluated  by  summing  over  the  particles  before  other  quantities  may 
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be  interpolated.  However,  Equation  3.23  allows  density  to  be  evolved  concurrently 
with  particle  velocities  and  other  field  quantities,  thus  significantly  reducing  the  com¬ 
putational  effort.  Although  Equation  3.23  does  not  conserve  mass  exactly  (Equation 
3.7  does,  provided  that  the  total  number  and  masses  of  particles  are  constant),  this 
is  usually  only  important  for  faster  flows  involving  shocks.  Direct  particle  summa¬ 
tion  can  be  used  intermittently  during  a  simulation  to  prevent  significant  "drift  in 
particle  densities  regarding  to  the  fluid  masses  around. 

3.2.3  Evaluation  of  Acceleration 

The  evaluation  of  particle  accelerations  comes  from  the  SPH  form  for  the  momen¬ 
tum  equation  which  is  dependent  on  how  to  evaluate  the  pressure  gradient  accelera¬ 
tion  and  implement  the  fluid  viscosity.  Benz  (1990)  had  detailed  derivation  for  this 
equation  for  inviscid  fluids. 

3. 2. 3.1  Pressure  Gradient  Acceleration 


In  SPH,  the  pressure  gradient  acceleration  term  in  Equation  3.15  is  usually  sym¬ 
metrized  by  writing, 


Vp 

P 


v  -  +  4vP. 


(3.26) 


This  results  in  the  most  common  SPH  expression  for  the  term, 

-(ivp)o  =  -£mt(|  +  g)v^  (3,7, 

Provided  that  the  kernel  is  an  even  function  of  r,  Equation  3.27  conserves  linear 
and  angular  momenta  exactly  since  the  forces  acting  between  individual  particles  are 
antisymmetric.  It  is  natural  because  of  the  law  of  action  and  reaction.  Momentum 
conservation  can  be  satisfied  by  an  infinite  number  of  symmetric  forms  of  the  pressure 
gradient  acceleration  term  given  by  Monaghan  (1992), 


-  Vp 


=  -^m6 

b 


Pa  Pb 


Pipl~e  PlPl~£ 


vgvtq6, 


(3.28) 
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where  s  may  take  any  value.  The  following  form  (i.e.,  e  —  1)  provides  ceitain  ad\an- 
tages  for  problems  involving  contact  discontinuities,  and  was  used  for  the  research 
described  herein. 

-(ivP)a=-£ny^)v^. 

3. 2. 3. 2  Viscosity 


Many  forms  of  artificial  viscosity  have  been  proposed  for  modeling  viscous  finid  in 
SPH  (Benz  1990:  Monaghan  1992).  The  most  commonly  used  form  is  incorporated 
into  the  momentum  equation  as. 


dua 

dt. 


b 


(  Pa+Pb 
\  Pa  Pb 


V0W'a6  +  g. 


where. 


{— 0.5a(cn  +  CbYfiab  +  3j£t 
0.5(pa  +  Pb) 

0 


if  ua6  •  Tab  <  0: 


otherwise, 


(3.30) 


(3.31) 


~  ab  ^  ab 

^  =  r2b  +  0.01  h2 


(3.32) 


This  viscosity  is  Galilean  invariant  and  conserves  total  linear  and  angular  momenta. 
The  na6  term  produces  a  shear  and  bulk  viscosity  which  permit  the  modeling  of  strong 
shocks,  a  and  0  are  chosen  to  be  1  and  2,  respectively,  for  best  results.  The  O.Ol/i2 


term  is  included  to  keep  the  denominator  nonzero.  In  this  form,  viscosity  is  intended 
only  to  provide  the  needed  dissipation  at  a  shock  front  to  convert  kinetic  energy  into 
internal  energy,  so  it  is  only  active  for  approaching  particles. 

Although  this  formulation  has  been  used  to  model  real  viscosity  (Artymowicz  and 
Lubow  1994),  it  produced  inaccurate  velocity  profiles  for  simulations  in  this  work. 
Equation  3.30  guarantees  conservation  of  angular  momentum,  which  is  important 
for  applications  involving  relatively  large  fluid  velocities  or  an  unbounded  fluid  edge. 
Since  the  applications  in  porous  media  flow  involve  low  velocities  and  SPH  particles 


fill  all  space,  a  more  realistic  form  of  viscosity  has  been  adopted.  Other  expressions 
have  been  proposed  to  model  real  viscous  forces:  however,  their  implementation  re¬ 
quires  nested  summations  over  the  particles  and  hence  twice  the  computational  effort 
(Flebbe  et  al.  1994;  Watkins  et  al.  1996).  or  second  derivatives  of  the  kernel  which 
introduce  sizable  errors  at  low  resolution  interpolations  (Takeda  et  al.  1994).  This 
method  employed  an  SPH  estimation  of  viscous  diffusion  as, 
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(^V"u)a  = 


Ef^biPa  d~  Pfc)rafe  ’  Tall  Qfe 

"  o  ■  ^  >  '■>'  Uq6' 


(3.33) 


PaPb{r2ab  +  0.01  h2) 

—  u 

Equation  3.33  is  based  on  a  similar  SPH  expression  used  by  Monaghan  (1995a)  to 
model  heat  conduction.  This  hybrid  expression  combines  a  standard  SPH  first  deriva¬ 
tive  with  a  finite  difference  approximation  of  a  first  derivative.  By  taking  a  Taylor 
expansion  about  particle  a,  it  can  be  shown  that  this  expression  is  approximate  (Mon¬ 
aghan  1995a).  This  formulation  conserves  linear  momentum  exactly,  while  angular 
momentum  is  only  approximately  conserved.  If  the  kernel  takes  the  form  of  Equation 
3.3.  Equation  3.33  simplifies  to, 

mb(na  +  Pb)uab  1  dWab 


(rV2u)0  =  y 
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(3.34) 


Substituting  Equations  3.29  and  3.34  into  Equation  3.15,  the  SPH  form  of  the  mo- 
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mentum  equation  is, 
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dt 
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V  PaPb  )  V  PaPt 


kail  drab 


+  g,  (3-35) 


which  is  used  to  evaluate  particle  accelerations. 


3.2.4  Dynamic  Pressure 

For  low  Reynolds  number  flows,  local  variations  in  pressure  gradient  which  force 
fluid  motion  can  be  very  small  in  comparison  to  the  hydrostatic  pressure  gradient. 
This  is  of  special  significance  to  SPH  since  pressure  is  obtained  using  an  explicit  func¬ 
tion  of  density  and  is  only  accurate  to  about  1%.  Consequently,  for  many  problems, 
it  is  simpler  to  model  the  dynamic  pressure  pa  defined  as, 


Pd  =  Pt~  Ph, 


(3.36) 


where  pt  and  ph  are  the  total  and  hydrostatic  pressures,  respectively.  Thus,  it  is 
the  dynamic  pressure  which  is  modeled  by  Equation  3.18  and  used  in  Equation  3.35. 
Since  pressure  appears  in  the  Navier-Stokes  equations  only  as  a  gradient,  the  effect 
of  ph  is  that  of  a  body  force. 

— -Vpt  =  — -Vp<f  —  -Vph-  (3.37) 

p  P  P 

Substituting  Equation  3.37  into  Equation  3.15  gives, 

—  =  --\/pd  +  //V2u  +  g - Vp*  —  — TJpd  +  i'V-u  +  F,  (3.38) 

dt  p  p  p 

where  F  is  the  net  body  force  driving  the  flow  defined  as. 

F  =  g  -  -Vph,  (3-39) 

P 

i.e..  fluid  particles  move  in  response  to  an  imbalance  of  forces  due  to  gravity  and  a 
large-scale  static  pressure  field.  Using  this  approach,  pressure  gradient  driven  flow 
through  a  periodic  lattice  can  be  easily  simulated.  For  simplicity,  p  is  used  in  the 
followings  to  denote  the  dynamic  pressure  Pd- 

3.2.5  Boundary  Conditions 

Initial  applications  of  quasi-incompressible  SPH  involved  high  Reynolds  number 
simulations  of  free  surface  flows  interacting  with  free-slip  boundaries  (Monaghan  1994; 
Monaghan  1995b;  Monaghan  1995c).  Such  work  employed  boundary  particles  which 
exerted  strong  repulsive  forces  to  prevent  SPH  particles  from  penetrating  solid  sur¬ 
faces.  The  force  has  the  form  of  Lennard-Jones  form  for  forces  between  molecules. 
To  realistically  model  porous  media  flow,  no-slip  fluid-solid  boundary  conditions  are 
needed.  In  addition,  for  the  free  surface  flows  considered  by  Monaghan  (1994),  bound¬ 
ary  particles  do  not  contribute  to  the  density  of  the  free  SPH  particles,  thus  permitting 
the  fluid  to  freely  leave  a  solid  boundary  with  no  pressure-driven  restoring  force.  In 
this  work,  boundary  particles  contribute  to  the  density  of  fluid  particles  such  that 
pressure  decreases  when  fluid  and  boundary  particles  diverge.  It  is  possible  to  imple¬ 
ment  such  a  boundary  condition  using  image  particles  (Libersky  and  Petschek  1990; 
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Libersky  et  al.  1993;  Randles  and  Libersky  1996).  The  images  are  created  by  re¬ 
flecting  fluid  particles  across  the  boundary  with  opposite  velocities.  This  procedure 
works  well  for  straight  channels,  but  introduces  density  errors  for  curved  surfaces. 
Takeda  et  al.  (1994)  achieved  a  no-slip  boundary  condition  using  special  boundary 
terms  which  mimic  a  half-space  filled  with  SPH  image  particles.  \\  hile  the  approach 
has  proved  useful  for  compressible  and  moderate  to  high  Reynolds  number  flows,  it 
did  not  yield  sufficiently  stable  results  for  quasi-incompressible  low  Reynolds  number 
simulations. 

In  this  work,  real  SPH  particles  are  used  to  create  no-slip  boundaries.  These 
particles  contribute  to  the  usual  SPH  expressions  for  density  and  momentum,  and 
their  own  densities  are  also  evolved.  Evolving  the  densities  of  boundary  particles  was 
found  to  better  capture  peak  pressures  than  if  the  densities  of  boundary  particles 
were  kept  constant.  Obstacles  are  created  by  initially  placing  SPH  particles  on  a 
regular  lattice  and  deleting  those  particles  which  fall  within  the  solid  matrixes.  Then 
boundary  particles  are  placed  in  layers  parallel  to  the  solid  surfaces.  Figure  3.2 
illustrates  the  concept  for  a  curved  boundary. 

For  each  fluid  particle  a,  the  normal  distance  da  to  the  boundary  is  calculated. 
This  normal  defines  a  tangent  plane  (a  line  in  two  dimensions)  from  which  the  normal 
distance  dB  to  each  boundary  particle  B  is  calculated.  The  velocity  of  particle  a  is 
extrapolated  across  the  tangent  plane,  assuming  zero  velocity  on  the  plane  itself, 
giving  each  boundary  particle  the  velocity, 

uB  =  ~~Tua-  (3-4°) 

da 

By  doing  so,  boundary  particles  are  assigned  artificial  velocities  such  that  antisymme¬ 
try  in  the  velocity  field  is  created  across  the  boundary  surface.  Ideally,  local  estimates 
of  the  velocity  gradients  at  the  surface  of  the  boundary  would  be  used  to  assign  these 
artificial  velocities  to  interior  points,  however,  such  estimates  would  require  a  second 
summation  over  the  particles  and,  hence,  a  substantial  increase  in  the  computation 
efforts.  The  approach  presented  here  is  simple,  stable,  accurate,  and  requires  little 
extra  computation. 


Boundary  particle 


Figure  3.2:  Construction  of  no-slip  fluid-solid  boundary  for  a  curved  surface. 

If  the  boundary  is  in  motion,  ua  should  be  replaced  by  the  fluid  velocity  rela¬ 
tive  to  the  boundary.  The  artificial  velocity  u#  is  used  to  calculate  viscous  foices, 
whereas  the  actual  boundary  velocity  is  used  to  evolve  boundary  particle  positions 
and  densities.  In  practice,  the  discrete  arrangement  of  boundary  particles  may  permit 
a  fluid  particle  to  closely  approach  the  nominal  curve  describing  the  boundary.  In 
such  circumstances,  potentially  large  artificial  velocities  for  boundary  particles  may 
result.  To  prevent  this  problem,  da  is  bounded  according  to. 

da  =  ma x(da,  -~^Ax),  (3.41) 

where  Ax  is  the  initial  nearest  neighbor  distance  between  fluid  particles  (Figure  3.3). 
Every  SPH  particle  has  its  mass  associated  with.  As  the  first  layer  of  boundary 
particles  is  placed  on  the  solid  surface,  the  porous  medium  simulated  actually  has  a 
lesser  porosity  than  intended.  This  effect  diminishes  as  Ax  decreases. 
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Fluid  volume  associated 
with  particle  a 


Figure  3.3:  Initial  hexagonal  configuration  of  SPH  fluid  particles. 


3.2.6  Interpolation  Kernel 


The  interpolation  kernel  is  used  to  calculate  a  weighted  sum  of  fluid  properties 
at  a  point.  The  use  of  different  kernels  for  SPH  is  analogous  to  the  use  of  different 
finite  difference  operators  for  the  finite  difference  method.  Most  SPH  applications 
employ  a  cubic  spline  kernel  (Schoenberg  1946;  Monaghan  and  Lattanzio  1985)  since 
it  resembles  a  Gaussian  while  having  compact  support, 


where, 


W(t,  h ) 


10 

77 rh2 


1  -  §s2  +  | s3 

<  (2-s)3 

1  3 


if  0  <  .9  <  1 ; 
if  1  <  s  <  2; 


0 


if  s  >  2, 


s  = 


h 


(3-42) 


(3.43) 


and  the  above  equation  is  normalized  for  two  dimensions.  However,  it  has  been 
shown  that  SPH  can  be  unstable  to  transverse  modes  when  kernels  with  compact 
support  are  used  (Morris  1994;  Morris  1996b).  As  higher-order  splines  more  closely 
approximating  a  Gaussian  are  employed,  these  instabilities  are  reduced.  One  reason 
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for  the  poor  performance  of  lower-order  splines  is  that  the  stability  properties  of  SPH 
depend  strongly  upon  the  second  derivative  of  the  kernel.  The  second  derivative  of  the 
cubic  spline  is  a  piecewise-linear  function,  and,  accordingly,  the  stability  properties  are 
inferior  to  those  of  smoother  kernels.  For  low  Reynolds  number  quasi-incompressiblo 
flow  simulations,  cubic  spline  kernel  rapidly  produced  significant  noise  in  pressure 
and  velocity  fields.  The  following  quintic  spline  kernel  (Schoenberg  1946)  has  been 
chosen  in  this  work,  here  normalized  for  two  dimensions,  which  results  in  much  less 
noise  in  the  pressure  and  velocity  fields. 


(3  _  sf  _  6(2  -  s)5  +  15(1  -  s)5  if  0  <  s  <  1: 


]V(r,h) 


i 
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I  (3  -  s)5  -  6(2  -  s)5 
1  (3  s)5 


0 


if  1  <  s  <  2; 
if  2  <  s  <  3; 
if  s  >  3. 


(3.44) 


Figure  3.4  shows  Wh2  as  a  function  of  s  for  the  quintic  spline  kernel  given  by  Equation 


Figure  3.4:  Wh?  versus  s  for  quintic  spline  kernel. 
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3.44.  By  using  a  quintic  spline,  particles  interact  within  a  distance  of  3 h  (Figure  3.1). 
This  produces  property  estimates  that  are  smoother  than  those  obtained  using  a 
cubic  spline.  The  numerical  solution  also  has  better  stability  properties  as  a  result 
(Morris  1996b).  The  quintic  spline  is.  however,  more  computationally  expensive  than 
the  cubic  spline  by  approximately  a  factor  of  two. 

3.2.7  Locating  the  Nearest  Neighboring  Particles 

As  defined  by  a  kernel  with  compact  support,  each  particle  has  a  finite  number 
of  "neighboring"  particles  with  which  it  interacts.  For  the  quintic  kernel,  particles 
interact  within  a  distance  of  3 h  (Figure  3.1).  The  problem  still  remains,  however,  to 
efficiently  identify  these  interacting  pairs  of  particles.  In  this  work,  an  algorithm  em¬ 
ploying  a  linked  list  data  structure  is  used  to  locate  neighboring  particles  (Monaghan 
and  Lattanzio  1985;  Morris  1996a). 

The  computational  domain  is  divided  into  square  cells,  each  with  a  side  length 
of  3 h.  A  list  of  particles  belonging  to  each  cell  is  created.  A  particle  located  within 
a  given  cell  then  must  only  consider  interactions  with  particles  in  neighboring  cells. 
The  lists  of  particles  within  each  cell  are  implemented  as  linked  lists.  That  is,  there 
is  a  pointer  to  the  first  particle  in  a  cell,  and  that  particle  then  points  to  the  second 
particle  and  so  on.  In  one  dimension,  the  algorithm  reads  as: 
for  i  =  1  to  all  particles  do 
j  =  int((xl  -  xmin) / (3h)) 
nexti  —  headj 
headj  =  i 

headj  is  a  pointer  to  the  first  particle  in  cell  j,  and  is  initially  set  to  0,  while  nexti  is 
a  pointer  from  particle  i  to  the  next  particle  in  the  linked  list,  headj  will  point  to  0 
if  cell  j  is  empty  and  nexti  will  point  to  0  if  particle  i  is  the  last  particle  in  the  list. 
Finding  the  nearest  neighbors  of  particle  i  in  cell  j  is  now  performed  as: 
for  cell  =  j  —  1  to  j  +  1  do 
k  —  headceii 
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while  (k  ^  0)  do 
consider  particle  k 
k  =  nextk 

However,  in  practice,  it  is  much  more  efficient  to  create  a  temporary  list  of  particles 
in  a  given  cell  and  interacting  neighboring  cells  and  evaluate  the  interactions  between 
them.  Also,  it  is  enough  to  only  consider  the  neighboring  cell  to  one  side  of  the  home 
cell  because  the  cell  on  the  other  side  will  consider  this  home  cell  as  its  neighboring 
cell  and  include  it.  For  example,  considering  the  SPH  interactions  for  cell  j.  a  list  ot 
particles  from  cell  j  (home  cell)  and  cell  j  +  1  (the  neighboring  cell  to  the  light)  is 
created: 

i  =  0 
k  =  headj 
while  ( k  ^  0)  do 
i  =  i  +  1 
listi  =  k 
k  =  nextk 

H'h.ome  —  ^ 
k  —  hc&dj+i 
while  ( k  ±  0)  do 
2=2  +  1 
listi  =  k 
k  =  nextk 

Tltotal  ^ 

Then  the  SPH  contributions  can  be  readily  obtained  by  considering  pairs  of  particles 

as: 

for  =  1  tO  Xifiome  tlo 

CL  =  ilStjj 

for  z2  =  H  to  ntotai  do 
b  =  listi2 

consider  particles  a  and  b 
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These  loops  consider  all  the  interactions  between  particles  in  the  home  cell  j  with 
each  other  and  with  those  particles  in  the  neighboring  cell  j  + 1.  The  first  loop  is  over 
rihome  particles  while  the  second  is  over  ntotai  particles  since  the  interactions  between 
particles  in  the  neighboring  cell  with  themselves  are  not  intended. 


3.2.8  Time  Integration 


The  time  integration  of  SPH  equations  can  be  performed  using  the  same  basic 
approaches  which  are  employed  for  other  explicit  hydrodynamic  methods.  The  chosen 
method  should  provide  high  order  accuracy  with  a  minimum  number  of  sweeps  over 
the  particles.  In  this  model,  explicit  time  integration  is  performed  using  the  modified 
Euler  method  in  which  the  time  step  r  is  limited  by  stability  constraints.  The  CFL 
condition  (Courant  et  al.  1928)  essentially  states  that-  the  maximum  numerical  rate 
of  propagation  of  information  must  exceed  the  physical  rate.  In  SPH.  this  translates 
to. 


t  <  0.25—.  (3.45) 

c 

Additional  constraints  arise  from  the  magnitude  of  particle  accelerations  fa  (Mon¬ 
aghan  1992). 


and  viscous  diffusions, 


r  <  0.25  min 

a 


(3.40) 


r  <  0.125—.  (3.47) 

v 

Equation  3.47  is  based  upon  the  usual  condition  for  an  explicit  finite  difference  method 
simulating  diffusion.  For  simulations  having  high  resolution  (small  h)  or  large  vis¬ 
cosity,  Equation  3.47  is  typically  the  dominant  time  constraint.  The  choice  of  kernel 
and  the  initial  arrangement  of  particles  influence  the  coefficients  in  Equations  3.45  to 
3.47.  In  particular,  different  splines  can  have  different  “effective”  resolution  lengths 
for  the  same  value  of  h.  For  example,  use  of  a  cubic  spline  (which  is  “narrower”  than 
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a  quintie  for  the  same  smoothing  length)  may  require  slightly  smaller  coefficients  in 
the  above  expressions. 

The  time  stepping  is  carried  out  using  a  predictor-corrector  scheme  (Al-Khafaji 
and  Tooley  1986).  The  following  equations  are  used  to  obtain  the  field  quantities  at 
the  next  time  step, 


^  = «°  +  ' 

(3.48) 

x3  =  x°  +  —  U°. 

(3.49) 

P*=P°  +  §# 

(3.50) 

=  J(m,x^,uKpK-), 

(3.51) 

I  _  I  __  1 

f  =  f(m:x2,u2, ...), 

(3.52) 

U2  =  U°  +  —  fa  , 

(3.53) 

I  n  7"  i 

X 2  =  X  +  -u2 , 

(3.54) 

(3.55) 

ul  =  2  —  u°, 

(3.56) 

x1  =  2x2  —  x°, 

(3.57) 

P 1  =  2p2  -  p°, 

(3.58) 

where  the  superscripts  are  the  time  step  index,  and  fa  and  fp  are  the  acceleration  and 

rate  of  change  of  density,  respectively.  In  practice,  the  predictor  step  uses  values  of  fa 

_  1  , 

and  fp  at  the  previous  midpoint,  e.g.,  /°  =  fa  2 ,  since  this  reduces  the  computation 
work  without  changing  the  order  of  the  scheme. 


84 


3.3  Model  Initialization  and  Execution 
3.3.1  Initialization 

For  each  simulation,  fluid  particles  are  initially  placed  throughout  the  compu¬ 
tational  domain,  i.e.,  the  unit  cell  of  the  periodic  porous  medium,  in  a  hexagonal 
arrangement  (Figure  3.3).  Each  particle  is  assigned  an  initial  density  p0.  an  initial 
velocity  of  zero,  and  a  mass  ma  =  po\ Q,  where  1 Q  is  the  volume  associated  with  parti¬ 
cle  a  (Figure  3.3).  Fluid  particles  that  fall  within  the  solid  matrixes  are  deleted.  The 
choice  for  the  number  of  fluid  particles  is  governed  by  the  desired  resolution  and  the 
computational  expense.  Larger  numbers  of  particles  produce  more  accurate  results 
but  also  increase  computation  time.  For  porous  media  simulations,  it  is  recommended 
that  each  pore  throat  should  be  spanned  by  at  least  15  particles.  As  a  result,  the 
computation  time  for  problems  having  small  porosity  (i.e.,  narrow  pore  throats)  dra¬ 
matically  increases  due  to  the  larger  number  of  particles  and  the  reduced  numerical 
time  step  (Equations  3.45  to  3.47). 

Boundary  particles  are  placed  on  the  solid  inclusions  using  a  pseuclo-hexagonal 
arrangement  (Figure  3.2).  The  first  layer  of  boundary  particles  are  positioned  on 
the  perimeter  with  an  equal  spacing  of  approximately  Ax  (depending  on  the  size  of 
the  solid  obstacles).  A  second  layer  of  an  equal  number  of  particles  is  then  placed 

/o 

at  a  distance  of  Ax  from  the  boundary.  A  third  layer  of  particles  is  positioned 
similarly  and  the  process  is  continued  until  boundary  particles  fill  an  annular  zone 
having  a  thickness  of  at  least  6h  from  the  solid  surface.  Particles  farther  than  6 h  are 
not  needed  because  they  would  not  contribute  to  the  calculations  due  to  the  use  of 
quintic  spline  kernel.  Once  the  boundary  particles  are  positioned,  each  is  assigned 
a  mass  consistent  with  its  contributing  volume  and  an  initial  density  of  po .  The 
solid  inclusion  is  usually  large  enough  to  accommodate  the  annular  zone  of  boundary 
particles.  However,  when  it  is  small,  Ax  may  be  reduced  to  enforce  the  boundary 
particle  condition  (i.e.,  an  annular  zone  having  a  thickness  of  at  least  6 h  from  the  solid 
surface).  In  this  work,  the  number  of  total  SPH  particles  Npart  (including  boundary 
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particles)  chosen  for  a  simulation  is  either  determined  by  having  about  20  particles 
spanning  the  narrowest  pore  throat  or  having  enough  boundary  particles  for  the  solid 
inclusions. 

The  positioning  of  boundary  particles  disrupts  the  initial  hexagonal  arrangement 
of  fluid  particles  in  the  problem  domain.  To  account  for  changes  in  particle  density 
near  the  boundaries,  the  initial  density  field  (fluid  particles  and  boundary  particles) 
is  recalculated  using  direct  summation  (Equation  3.7).  Once  particle  densities  are 
corrected,  a  body  force  F  is  applied  to  initiate  fluid  motion.  F  is  related  to  the 
hydraulic  gradient  i  by. 

i  =  -.  (3.59) 

9 

In  this  model,  F  can  be  applied  in  any  direction.  Periodic  unit  cell  boundary  condi¬ 
tions  are  applied  to  all  fluid  quantities  and  the  flow  is  driven  by  the  effective  body 
force. 

3.3.2  Execution 

Once  the  body  force  is  applied  to  the  system,  particle  densities  are  evolved  accord¬ 
ing  to  Equation  3.23  and  particle  accelerations  are  computed  using  Equation  3.35. 
The  smoothing  length  h  is  chosen  equal  to  2 Ax  for  the  quintic  kernel.  To  ensuie  sta¬ 
bility  of  the  integration  scheme,  the  time  step  is  limited  according  to  the  conditions  set 
forth  in  section  3.2.8.  Particle  velocities,  positions,  and  densities  are  updated  using 
a  predictor-corrector  method  (section  3.2.8).  No-slip  fluid-solid  boundary  conditions 
are  simulated  by  assigning  artificial  velocities  to  the  boundary  particles  using  Equa¬ 
tion  3.40.  Although  computationally  less  expensive,  Equation  3.23  does  not  conserve 
mass  exactly  and  may  introduce  errors  over  the  course  of  a  simulation.  To  correct  foi 
this',  particle  densities  are  updated  every  50  time  steps  by  direct  summation  (Equa¬ 
tion  3.7).  Using  this  procedure,  computational  speed  is  maintained  and  fluid  mass  is 
conserved  as  well.  Locating  neighboring  particles  is  achieved  by  creating  the  linked 
lists  for  SPH  particles  (section  3.2.7). 


86 


Periodic  unit  cell  boundaries  are  created  through  the  use  of  solid  inclusion  and 
fluid  particle  images  (Figure  3.5),  and  by  wrapping  fluid  particles  around  the  flow 
domain  when  they  leave  one  boundary  (Figure  3.6).  Image  fluid  particles  and  solid 
inclusions  are  created  within  a  distance  of  3/z  (i.e..  one  quintic  kernel  radius)  from 
the  unit  cell  to  provide  the  necessary  "neighbors"  for  particles  within  the  cell. 


Boundary  for  image  particles 
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Unit  cell  boundary 


Fluid  particle 


Figure  3.6:  Wrapping  fluid  particles  around  unit  cell. 


respectively.  In  SPH,  they  are  evaluated  as, 


v  unit  cell  ^  Pb 


Vy  Vunit  cell  ^ 


(3.63) 


where  uj>x  and  are  x  and  y  components  of  the  velocity  of  particle  b,  lespectivelv. 
During  the  course  of  a  simulation,  vx,  vy,  and  the  maximum  particle  velocity  umax 
are  recorded.  Changes  in  these  values  with  time  are  used  to  determine  whether  or 
not  the  flow  has  reached  steady  state. 


3.4  Flow  Model  Verification 

The  SPH  numerical  flow  model  has  been  tested  in  many  cases.  Simulations  using 
the  method  show  close  agreement  with  series  solutions  for  Couette  and  Poiseuille 
flows  and  with  other  solutions  for  flow  past  regular  lattices  of  obstacles. 
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3.4.1  Simulations  of  Couette  and  Poiseuille  Flows 

3. 4. 1.1  Couette  Flow 


Couette  flow  is  simulated  between  infinite  plates  located  at  y  =  0  and  y  =  L 
(Figure  3.7).  The  system  is  initially  at  rest.  At  time  /  =  0.  the  upper  plate  moves  at 


constant  velocity  uo  parallel  to  the  x  axis.  The  series  solution  for  the  time-dependent 
behavior  of  this  flow  is, 


ux(y,t) 


u°  ,  2U°  /  . 

—  y  +  >  - (  — 1)  sin 

L  '  mr 

n=  1 


(3.64) 


where  ux  is  the  fluid  velocity  in  the  x  direction.  The  flow  was  simulated  using  SPH 

2  h 

for  v  =  10~6  ^5-.  L  =  10_3m,  p  =  103  — tin  =  1.25  x  10~3  ^r.  and  with  50  particles 
*  '  m  * 

spanning  the  channel.  This  corresponds  to  a  Reynolds  number  of  0.0125  if  using, 


Re  =  ^.  (3.65) 

u 

Figure  3.8  shows  a  comparison  between  velocity  profiles  obtained  using  Equation  3.64 
and  SPH  at  several  times  including  the  steady  state  solution  ( t  —  oc).  The  results  are 
in  close  agreement  within  0.5%,  confirming  the  accuracy  of  the  approach  used  to  eval¬ 
uate  viscous  and  boundary  forces  with  SPH.  Lower  resolution  simulation  completed 
with  20  particles  spanning  the  channel  was  found  to  agree  to  within  approximately 
2%  of  the  series  solution  values. 


y  (1 0”3  m) 


Figure  3.8:  Comparison  of  SPH  and  series  solutions  for  Couette  flow  (Re  -  0.0125). 


3.4. 1.2  Poiseuille  Flow 

Poiseuille  flow  between  stationary  infinite  plates  at  y  =  0  and  y  =  L  is  simulated 
(Figure  3.9).  The  fluid  is  initially  at  rest  and  driven  by  an  applied  body  force  F 


Figure  3.9:  Poiseuille  flow. 
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Figure  3.10:  Comparison  of  SPH  and  series  solutions  for  Poiseuille  flow  (Re.  =  0.0125). 


parallel  to  the  x  axis  for  t  >  0.  The  series  solution  for  the  transient  behavior  is. 


ux(y.  t ) 


F_ 

2u 


y(L  -  y)  +  Y_ l 


AFL2 


z/7r3(2n  +  l)3 


sin 


(? 


(2 n  +  1))  exp  ^ 


(2 n  +  l)2-^2!/  \ 

Z5  V‘ 

(3.66) 


2  7, 

SPH  was  used  to  simulate  Poiseuille  flow  for  u  =  10-6  L  =  10-3  m.  p  =  103  — 

s  m3 

F  =  10~4  and  with  50  particles  spanning  the  channel.  This  corresponds  to  a 
s 

peak  fluid  velocity  it0  =  1.25  x  10-5  ^  and  a  Reynolds  number  of  0.0125  if  using 
Equation  3.65.  A  comparison  between  velocity  profile  obtained  using  Equation  3.66 
and  SPH  appears  in  Figure  3.10.  The  results  are  again  in  close  agreement  with  the 
largest  discrepancy  being  about  0.7%  for  the  steady  state  solution.  Lower  resolution 
simulation  completed  with  20  particles  spanning  the  channel  was  found  to  agree  to 
within  approximately  2%  of  the  series  solution  values. 


3.4.2  Simulations  of  Flow  Through  Periodic  Lattices  of  Obstacles 

The  Couette  and  Poiseuille  flow  simulations  tested  the  interaction  between  viscous 
and  body  forces  and  the  effectiveness  of  the  no-slip  boundary  conditions  in  the  model. 
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However,  these  flows  are  essentially  one-dimensional  and  do  not  produce  variations 
in  dynamic  pressure.  A  more  challenging  test  of  the  method  involves  flow  through 
lattices  of  obstacles.  In  this  SPH  model,  circular  and  elliptical  obstacles  ha\o  boon 
implemented.  Before  presenting  the  simulations,  three  kinds  of  regular  arrays  of 
obstacles  for  spatially  periodic  porous  media  considered  throughout  this  report  are 
introduced.  These  are  square,  staggered,  and  hexagonal  arrays  of  obstacles. 

3.4. 2.1  Square,  Staggered,  and  Hexagonal  Arrays  of  Obstacles 

According  to  Perrins  et  al.  (1979),  there  are  five  possible  ways  of  packing  obsta¬ 
cles  in  regular  two-dimensional  arrays  (International  tables  for  X-ray  crystallography 
1962).  Three  of  these  arrays — square,  staggered,  and  hexagonal — are  used  in  this 
work.  Figures  3.11  to  3.13  depict  the  unit  cell  geometry  for  the  three  arrays  with 
circular  and  elliptical  inclusions.  The  unit  cell  for  the  square  or  staggered  array 


L 


~\ 


L 


d 


Figure  3.11:  Two-dimensional  cross  sectional  representations  of  unit  cell  geometry  for 
square  array  with  circular  and  elliptical  inclusions. 
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Unit  cell  boundary 


L  i  .  L 


Figure  3.12:  Two-dimensional  cross  sectional  representations  of  unit  cell  geometry  for 
staggered  array  with  circular  and  elliptical  inclusions. 

has  side  length  L,  and  the  unit  cell  for  the  hexagonal  array  has  side  lengths  L\  and 
L>  =  x/3Li.  The  staggered  array  may  be  regarded  as  a  superposition  of  two  square 
arrays.  The  circular  inclusion  is  simply  characterized  by  the  radius  R,  while  the  ellip¬ 
tical  inclusion  is  characterized  by  the  major  axis  2a,  minor  axis  2b.  and  the  angle  a 
(0°  <  a  <  180°)  between  the  direction  of  the  major  axis  and  the  positive  x  direction. 
For  circular  inclusions,  the  three  lattices  are  isotropic  for  fields  applied  in  the  plane, 
while  for  elliptical  inclusions,  intrinsic  anisotropy  exists. 

3. 4. 2. 2  Flow  Through  a  Periodic  Square  Lattice  of  Circular  Cylinders 

The  flow  through  a  periodic  square  lattice  of  circular  cylinders  was  first  simulated. 
This  particular  configuration  has  been  studied  extensively  as  a  simple  model  of  flow 
through  fibrous  porous  media.  The  periodic  flow  was  simulated  using  SPH  for  L  = 
0.1  m,v  =  10-6  R  =  0.02m,  F  =  1.5 x  10-7  and  c  =  5.77x  10-4  Replacing 
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L  with  R  in  Equation  3.65  and  taking  the  velocity  scale  to  be  it0  =  5  x  10“°  ™ 
gives  Re  =  1.  The  SPH  simulation  was  run  using  approximately  3000  particles 
with  a  nearest  neighbor  distance  Ax  of  0.002  rn.  Steady  state  was  reached  after 
approximately  1500  time  steps.  To  investigate  long-term  behavior,  the  simulation 
was  continued  for  another  6000  time  steps  such  that  particle  arrangements  became 
disordered.  The  problem  was  also  modeled  using  a  finite  element  method  (FE.M) 
program  for  steady  incompressible  viscous  flow.  Velocity  and  pressure  distributions 
from  the  two  solutions  were  compared  by  plotting  values  within  one  nearest  neighbor 
distance  of  the  four  paths  described  in  Figure  3.14.  The  results  were  also  compared 


Unit  cell 


Figure  3.14:  Paths  for  comparison  of  SPH  and  FEM  solutions  for  flow  through  a 
periodic  square  lattice  of  circular  cylinders. 

using  contour  plots.  As  the  FEM  employs  a  mesh  to  obtain  a  solution,  it  is  relatively 
easy  to  obtain  contour  plots.  The  corresponding  plots  for  SPH  were  obtained  by 
interpolating  the  particle  quantities  to  a  50  by  50  array  of  grid  points  using  the 
quintic  kernel.  Smoothing  lengths  of  1  and  3  grid  spacings  were  used  for  the  velocity 
and  pressure,  respectively.  A  greater  amount  of  smoothing  was  needed  to  remove 
small  fluctuations  from  the  pressure  field.  Contours  generated  using  this  method  are 
inaccurate  in  the  immediate  vicinity  of  the  cylinder. 

Figure  3.15  shows  a  comparison  of  velocity  profiles  obtained  using  SPH  and  FEM 
for  paths  1  and  2  defined  in  Figure  3.14.  The  results  obtained  using  SPH  are  in  close 
agreement  with  those  from  FEM  throughout  the  flow  domain.  Corresponding  contour 
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plots  of  velocity  magnitude  are  shown  in  Figure  3.16.  Good  agreement  is  obtained  for 
the  bulk  of  the  flow,  although  the  contour  smoothing  method  inaccurately  represents 
SPH  velocities  near  the  cylinder. 

Figure  3.17  shows  the  dynamic  pressure  along  paths  3  and  4  in  Figure  3.14.  The 
arc  of  path  3  was  taken  as  0.002  m  (one  SPH  nearest  neighbor  distance)  beyond  the 
cylinder  boundary.  The  SPH  dynamic  pressure  profile  shows  small  local  fluctuations 
in  the  vicinitv  of  the  cylinder.  Elsewhere,  the  SPH  and  FEM  solutions  are  in  close 
agreement.  The  peaks  in  the  pressure  obtained  using  SPH  on  the  boundaic  itself  fell 
short  of  the  FEM  results  by  approximately  ST.  The  FEM  better  captures  pressure 
extrema  since  grid-stretching  increases  resolution  in  the  vicinity  of  the  cylinder.  Cor¬ 
responding  contour  plots  of  pressure  given  in  Figure  3.18  show  that  good  agreement 
is  again  obtained  for  the  bulk  of  the  flow.  This  simulation  was  repeated  with  twice 
the  particle  resolution  (approximately  11000  particles)  and  peak  pressures  were  re¬ 
produced  to  within  5%.  Pressure  values  in  the  immediate  vicinity  of  the  boundary, 
however,  still  exhibited  small  fluctuations. 

Flow  through  a  periodic  square  lattice  of  circular  cylinders  was  also  solved  for 
L  =  0.1m,  v  =  lO-4^,  R  =  0.02m,  F  =  5  x  10“5  and  c  =  1  x  10“2 
Uo  =  1.5  x  10“4  ™  gives  Re  =  0.03.  The  steady  state  of  this  flow  was  reached  after 
approximately  1500  time  steps.  However,  the  initial  lattice  was  relatively  unchanged 
at  this  point.  To  demonstrate  the  long-term  behavior  of  the  method,  the  simulation 
was  run  for  another  300000  time  steps.  A  comparison  of  velocities  along  paths  1 
and  2  appears  in  Figure  3.19  and  velocity  contour  plots  are  shown  in  Figure  3.20. 
The  results  obtained  using  SPH  are  in  close  agreement  with  those  obtained  by  FEM. 
A  comparison  oCpressure  fields  presented  in  Figures  3.21  and  3.22  also  shows  close 
agreement  for  the  bulk  of  the  flow,  with  similar  fluctuations  as  observed  for  Re  =  1. 
Although  pressure  extrema  on  the  boundary  were  not  fully  captured  by  SPH,  the 
discrepancies  were  somewhat  smaller  (about  5%).  The  simulation  was  repeated  foi 
fewer  time  steps  with  twice  the  resolution  (approximately  11000  particles)  and  peak 
pressures  were  obtained  with  less  than  4%  discrepancy.  Once  again,  however,  small 
pressure  fluctuations  were  observed  near  the  boundary. 


Figure  3.15:  Comparison  of  SPH  and  FEM  velocity  profiles  along  paths  1  and  2  in 
Figure  3.14  for  Re  —  1. 
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Figure  3.16:  Contour  plots  of  velocity  magnitude  using  (a)  FEM,  and  (b)  SPH  for 
Re  =  1  (contour  lines  are  labeled  in  unit  of  10-4  ^r). 


Figure  3.17:  Comparison  of  SPH  and  FEM  pressure  profiles  along  paths  3  and  4  in 
Figure  3.14  for  Re  =  1. 
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Figure  3.18:  Contour  plots  of  pressure  using  (a)  FEM,  and  (b)  SPH  for  Re  = 
(contour  lines  are  labeled  in  unit  of  10-6  Pa). 
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Figure  3.19:  Comparison  of  SPH  and  FEM  velocity  profiles  along  paths  1  and  2  in 
Figure  3.14  for  Re  =  0.03. 
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Figure  3.20:  Contour  plots  of  velocity  magnitude  using  (a)  FEM,  and  (b)  SPH  for 
Re  —  0.03  (contour  lines  are  labeled  in  unit  of  10-4  ^r). 
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Figure  3.22:  Contour  plots  of  pressure  using  (a)  FEM,  and  (b)  SPH  for  Re 
(contour  lines  are  labeled  in  unit  of  10~3  Pa). 


an 
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3. 4. 2. 3  Flow  Through  a  Periodic  Square  Lattice  of  Elliptical  Cylinders 

The  flow  through  a  periodic  square  lattice  of  elliptical  cylinders  was  also  simulated 
using  both  SPH  and  FEM  for  comparative  study.  Velocity  and  pressure  distributions 
from  the  two  solutions  were  compared  by  plotting  values  within  one  nearest  neighbor 
distance  of  the  four  paths  described  in  Figure  3.23.  The  problem  was  solved  for 


Unit  cell 


Figpre  3.23:  Paths  for  comparison  of  SPH  and  FEM  solutions  for  flow  through  a 
periodic  square  lattice  of  elliptical  cylinders. 

L  =  2.16mm,  v  =  lO-6^^,  a  =  0.771  mm,  Q.  =  2,  a  =  45°,  F  =  0.001  and 

5  b  s 

c  =  0.012  ^r.  The  steady  state  Darcy  velocity  vx  was  found  to  be  8.10  x  10“°  ^  using 
SPH  and  8.29  x  1CT5  ^  using  FEM.  The  two  values  agree  within  a  2.3%  difference. 
The  average  value  of  vx  gives  Re  =  0.06  if  a  is  used  as  the  length  scale.  The  SPH 
simulation  was  run  for  50200  time  steps  using  41412  particles  with  a  nearest  neighbor 
distance  Ax  of  0.0108  mm. 

Figure  3.24  shows  a  comparison  of  velocity  profiles  obtained  using  SPH  and  FEM 
for  paths  1  and_2.  The  results  obtained  using  SPH  are  in  close  agreement  with  those 
from  FEM  throughout  the  flow  domain.  Figure  3.25  shows  a  comparison  of  dynamic 
pressure  along  paths  3  and  4.  The  two  solutions  are  in  close  agreement  except  at  the 
peaks  in  the  pressure  in  the  vicinity  of  the  cylinder.  The  fluctuations  are  similar  to 
those  observed  for  circular  inclusions  (Figures  3.17  and  3.21).  Figures  3.24  and  3.25 
confirm  that  SPH  is  capable  of  modeling  flow  past  elliptical  inclusions  in  general. 
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3. 4. 2.4  Flow  Through  Periodic  Square  and  Hexagonal  Lattices  of  Circular 
Cylinders 


The  SPH  model  was  further  used  to  simulate  two-dimensional  creeping  Hows  of 

i  2 

water  (p  =  103  =  10-6  --)  through  two  spatiallv  periodic  porous  media.  Tilt' 

m  b 

primary  objective  of  the  simulations  was  to  verify  the  results  of  the  model  through  a 
comparison  with  other  available  solutions  in  the  literature.  The  media  are  composed 
of  circular  cylinders  having  radius  R  that  are  arranged  in  regular  square  and  hexagonal 
lattices  (Figures  3.11  and  3.13).  For  the  square  lattice,  simulations  were  performed 
for  L  —  1.2  mm  and  porosity  n  =  0.3,  0.4.  0.5.  and  0.6.  The  corresponding  values  of 
cylinder  radius  are  R  =  0.566,  0.524,  0.479,  and  0.428  mm.  respectively.  Simulations 
for  the  hexagonal  lattice  were  performed  for  L\  —  1.2  mm  and  n  —  0.3.  0.4.  0.5.  and 
0.6.  The  corresponding  values  of  cylinder  radius  for  these  simulations  are  R  =  0.527. 
0.4S8.  0.446.  and  0.399 mm,  respectively.  The  cylinder  diameters  used  in  this  study 
fall  in  the  size  range  of  medium  sand.  In  particular,  to  have  a  better  view  about  the 
SPH  fluid  particles,  the  unit  cell  chosen  here  for  the  square  lattice  is  different  from  the 
one  in  Figure  3.11.  The  unit  cells  used  in  this  section  and  the  paths  for  comparison 
are  shown  in  Figure  3.26. 

During  the  course  of  a  simulation,  the  cylinder  drag  force  in  the  :r  direction,  i.e., 
the  direction  of  the  applied  body  force,  was  also  recorded  in  addition  to  vT  and  uxjnax 
every  100  time  steps.  The  cylinder  drag  force  in  the  x  direction  Fd  is  the  sum  of 
boundary  particle  forces  upon  the  cylinder  in  that  direction.  A  dimensionless  drag 
force  Fd  is  defined  as, 


Fd  = 


Fd_ 
pvx ' 


The  Reynolds  number  and  Mach  number  of  a  flow  are, 

2  Rvx 
Re  ~ - , 


and, 


(3.67) 


(3.68) 


Ma  = 


Vx_ 

7 

c 


(3.69) 
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respectively.  Tables  3.1  and  3.2  present  a  summary  of  the  results  obtained  for  the 

steady  state  flows.  The  small  values  of  Re.  correspond  to  creeping  flow  conditions  and 

the  small  values  of  Ma  indicate  that  each  flow  is  quasi-incompressible. 

For  any  SPH  simulation,  the  fluid  particles  will  eventually  become  disordered 

as  they  move  through  the  flow  domain.  Due  to  the  low  Reynolds  numbers,  however. 

steady  state  conditions  were  achieved  in  a  considerably  shorter  time  for  the  numerical 

simulations  in  this  study.  The  long-term  performance  of  the  model  was  investigated 

for  each  cvlincler  lattice  for  n  =  0.5  and  F  =  0.049  4i.  The  initial  and  final  particle 

s 

positions  are  shown  in  Figures  3.27  and  3.28  for  the  square  and  hexagonal  lattices, 
respectively.  The  model  was  run  for  39000  time  steps  for  the  square  lattice  and  27506 
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Figure  3.27:  Particle  positions  for  square  lattice  with  n  =  0.5  and  F  =  0.049  ^  (a) 
initial  positions^  and  (b)  final  positions.  Fluid  and  boundary  particles  are  shown  in 
black  and  grey,  respectively. 

steps  for  the  hexagonal  lattice.  In  each  case,  the  fluid  particles  were  fully  disordered 
at  the  end  of  the  simulation.  Figure  3.29  shows  uXyTnax  and  Fd  as  a  function  of  time 
for  each  simulation.  The  values  quickly  reached  steady  state,  and  thereafter  exhibited 
only  minor  fluctuations.  The  corresponding  plot  of  Darcy  velocity  is  shown  in  Figure 


(a)  (b) 
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Figure  3.28:  Particle  positions  for  hexagonal  lattice  with  n  =  0.5  and  F  =  0.049  -p 
(a)  initial  positions,  and  (b)  final  positions.  Fluid  and  boundary  pai tides  aie  shown 


in  black  and  grey,  respectively. 
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Figure  3.29:  Maximum  fluid  velocity  and  dimensionless  cylinder  drag  force  versus 
time  for  square  and  hexagonal  lattices  ( n  =  0.5,  F  =  0.049 

s 

3.30.  Steady  state  Darcy  velocity  was  also  achieved  quickly.  Figures  3.29  and  3.30 
illustrate  that  the  model  is  numerically  stable  well  after  steady  state  conditions  are 
reached.  Each  problem  was  also  solved  using  FEM  for  stead}-  incompressible  viscous 
flow.  In  terms  of  computation  time,  SPH  reached  steady  state  faster  than  FEM  for 
a  similar  number  of  particles  and  nodes.  However,  a  much  longer  time  was  required 
for  the  SPH  fluid  particles  to  become  fully  disordered.  Velocity  and  dynamic  pres¬ 
sure  distributions  from  the  two  methods  are  compared  by  plotting  values  within  one 
nearest  neighbor  distance  (Ax)  of  the  four  paths  (A,B,C,D)  shown  in  Figure  3.26. 
Figure  3.31  shows  profiles  of  fluid  velocity  in  the  x  direction  and  dynamic  pressure 
for  paths  A  and^B,  respectively,  through  the  square  lattice  unit  cell  for  n  =  0.5  and 

F  =  0.049  .  Corresponding  plots  for  paths  C  and  D  through  the  hexagonal  lattice 

s 

unit  cell  are  shown  in  Figure  3.32.  For  both  lattices,  the  values  are  in  close  agree¬ 
ment  with  typical  errors  of  about  5%.  The  plots  show  that  the  SPH  results  are  less 
smooth  than  those  obtained  using  the  FEM.  This  “noise"  is  due  to  particles  moving 
past  each  other,  which  is  then  amplified  by  the  relatively  stiff  equation  of  state.  In 
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Figure  3.30:  Darcy  velocity  versus  time  for  square  and  hexagonal  lattices  (/?  - 

0.5.  F  =  0.049  ^k). 

s 

general,  values  obtained  using  SPH  show  more  variability  as  the  compressibility  of 
the  fluid  decreases.  The  plots  also  indicate  that  computed  pressures  for  the  hexago¬ 
nal  lattice  are  not  fully  realized  near  the  stagnation  points  of  the  cylinders.  This  is 
consistent  with  the  findings  in  the  last  few  sections  which  showed  that,  at  comparable 
resolutions,  SPH  does  not  fully  capture  pressure  extrema  on  solid  boundaries.  The 
FEM  better  captures  the  extrema  since  grid-stretching  increases  resolution  in  the 
vicinity  of  a  cylinder.  The  SPH  and  FEM  solutions  are,  however,  in  close  agreement 
for  the  bulk  of  the  flow. 

Four  additional  simulations  were  performed  for  each  cylinder  lattice  for  n  =  0.5 
and  F  =  0.0392,  0.0245,  0.0098,  and  0.0049  Figure  3.33  shows  steady  state  Darcy 
velocity  varies  linearly  with  hydraulic  gradient  for  each  lattice,  which  is  in  agreement 
with  Darcy’s  law.  The  computed  values  of  hydraulic  conductivity  are  kHx  =  0.0255  ™ 
for  the  square  lattice  and  knx  =  0.0307  for  the  hexagonal  lattice.  At  n  =  0.5,  the 
hexagonal  lattice  has  a  higher  hydraulic  conductivity  than  the  square  lattice  because 
of  the  larger  flow  channels  in  the  hexagonal  lattice  unit  cell.  Values  of  kHx  are  plotted 
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Figure  3.31:  Comparison  of  SPH  and  FEM  results  for  square  lattice  (a)  velocity  profile 
for  path  A,  and  (b)  dynamic  pressure  profile  for  path  B  (n  =  0.5,  F  =  0.049  p1). 
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Figure  3.32:  Comparison  of  SPH  and  FEM  results  for  hexagonal  lattice  (a)  velocity 
profile  for  path  C,  and  (b)  dynamic  pressure  profile  for  path  D  (n  =  0.5,  F  =  0.049  3). 
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Figure  3.33:  Darcy  velocity  versus  hydraulic  gradient  for  square  and  hexagonal  lattices 
(n  =  0.5). 


as  a  function  of  n  for  both  lattices  in  Figure  3.34.  Noting  the  linearity  of  plots  in 
Figure  3.33.  values  of  kHx  for  media  having  porosities  other  than  0.5  were  calculated 
using  F  =  0.049  As  n  increases,  knT  increases  for  both  lattices  and  the  difference 
in  kHx  for  the  lattices  decreases. 

Values  of  and  Fj  for  each  numerical  simulation  are  shown  as  functions  of 
R 

porosity  in  Figure  3.35.  Solutions  obtained  using  FEM  as  well  as  published  results 

from  previous  studies  (Sangani  and  Acrivos  1982b;  Sangani  and  Yao  1988b;  Meegoda 

et  al.  1989;  Edwards  et  al.  1990;  Ghaddar  1995),  which  were  obtained  for  various 

values  of  R ,  are  also  shown  for  comparison.  One  additional  solution  has  been  included 

in  Figure  3.35(a)  for  the  square  lattice  at  n  =  0.558  for  comparison  with  the  data 

of  Meegoda  et  al.  (1989).  Values  are  in  good  agreement  with  a  maximum  difference 

of  57c,  indicating  that  SPH  is  capable  of  producing  results  that  are  comparable  with 

k 

those  obtained  using  other  methods.  Figure  3.35  also  indicates  that  is  a  unique 
function  of  n  for  each  periodic  lattice. 

Recalling  the  definition  of  in  Equation  3.67,  Equation  2.106  leads  to  the  fol¬ 
lowing  relation  for  the  media  considered  in  this  work, 
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Figure  3.34:  Hydraulic  conductivity  versus  porosity  for  square  and  hexagonal  lattices. 


‘^unit  cell  '  n5/cxFj.  (3.  /  0) 

where  ns  is  the  number  of  solid  obstacles  in  the  unit  cell  and  Aunitceii  is  the  area  of  the 
unit  cell.  This  relationship  provides  another  independent  check  of  SPH  simulation 
results,  which  appears  in  Tables  3.3  and  3.4.  It  is  evident  that  the  maximum  error  of 
the  SPH  numerical  results  is  about  6%. 

3.5  Summary 

Necessary  extensions  have  been  implemented  and  tested  which  allow  SPH  to  model 
incompressible  flow  through  porous  media.  Test  results  confirm  that  the  proposed 
modifications  to_the  equation  of  state,  viscosity  formulation,  boundary  conditions, 
and  interpolation  kernel  result  in  a  method  which  is  stable  and  accurate.  Simulations 
of  flow  through  spatially  periodic  porous  media  show  Darcy  velocity  proportional 
to  hydraulic  gradient,  as  required  by  Darcy’s  law.  In  addition,  the  solutions  are  in 
close  agreement  with  values  obtained  using  the  finite  element  method  and  published 
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Figure  3.35:  -=£  and  dimensionless  cylinder  drag  force  versus  porosity  for  square  and 
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hexagonal  lattices  (a)  -=&,  and  (b)  dimensionless  cylinder  drag  force. 
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CHAPTER  4  DEVELOPMENT,  VERIFICATION,  AND 
APPLICATION  OF  A  PORE-SCALE  DIFFUSION  MODEL  USING 
SMOOTHED  PARTICLE  HYDRODYNAMICS 

In  this  chapter,  the  method  of  SPH  is  formulated  to  solve  the  convection-diffusion 
equation  to  model  tracer  diffusion  through  porous  media.  Solutions  obtained  using 
SPH  are  compared  with  other  solutions  and  the  model  is  used  to  calculate  diffusion 
coefficients  of  spatially  periodic  porous  media  for  the  steady  state  diffusion  problem. 
Diffusion  coefficients  are  then  used  to  calculate  nondimensional  diffusivities  of  the 
media.  The  effects  of  media  properties  on  the  values  of  nondimensional  diffusivitv 
are  also  studied. 


4.1  Model  Establishment  for  Pore-Scale  Diffusion  Problem 


For  the  tracer  diffusion  problem,  concentration  C  becomes  a  property  associated 
with  SPH  particles  in  addition  to  mass  m  and  density  p.  If  the  convection-diffusion 
equation  (Equation  2.57)  is  rewritten  as, 

(4.1) 


—  =  d0V2C  =  ^V2C, 
dt  p 


it. 


becomes  evident  t licit  the  method  used  to  treat  viscosity  (Equation  3.33)  can  be 


adopted  to  evolve  tracer  concentration, 


dCa 

dt 


( d$p  do,bPb)Tab  ’  ab  ^  ^  ^  2) 

-\P  )a  h  PaPbir\b  ~b  0.01/i2) 


where. 


Cab  —  Ca  —  Cb- 


(4.3) 


Accordingly,  the  time  step  for  the  integration  scheme  is  limited  by, 


h 2 

t  <  0.125—. 

do 


(4.4) 


The  following  predictor-corrector  equations  are  used  to  obtain  particle  concentrations 
at  the  next  time  step, 

<^=C°  +  ^/°,  (4.5) 

fc  =  (T6) 

C 2  =  C°  -I-  -  ff.  ,  ( 4 .  i ) 

C1=2C5-C°,  (4.S) 

where  fc  is  the  rate  of  change  of  concentration  in  time.  As  in  the  flow  problem,  the 

_  i 

predictor  step  uses  the  value  of  fc  at  the  previous  midpoint  (i.e.,  /c  =  fc')- 

The  impenetrable  solid  boundary  condition  (Equation  2.131)  locally  requires  = 
0,  which  could  be  achieved  by  assigning  an  appropriate  artificial  concentration  for 
each  boundary  particle.  However,  this  approach  may  not  prevent  solute  mass  from 
diffusing  into  the  solid  phase  in  a  global  sense.  In  this'  model,  the  impenetrable  solid 
boundary  condition  is  easily  implemented  by  allowing  the  exchange  of  solute  mass 
only  between  fluid  particles. 

4.2  Calculation  of  Diffusion  Coefficients  of  Porous  Media  Using  SPH 

A  steady  state  diffusion  problem  is  solved  using  SPH  to  illustrate  the  calculation 
of  diffusion  coefficients  for  spatially  periodic  porous  media  (Figure  4.1).  Faces  x  = 
x\  and  x  =  x2  of  the  medium  are  maintained  at  constant  concentrations  C\  and 
C2,  respectively,  with  C\  >  Ci-  As  a  consequence  of  the  concentration  difference 
C\  —  C2,  a  steady  state  flow  of  solute  mass  necessarily  occurs  across  the  medium.  If 
the  molecular  diffusion  coefficient  do  is  assumed  to  be  constant,  that  is,  independent 
of  local  concentration  C  and  position  R,  it  is  expected  on  a  physical  basis  that  the 
concentration,  averaged  at  a  length  scale  greater  than  the  size  of  the  unit  cell,  will 
decrease  linearly  across  the  medium.  This  linear  behavior  would  be  obtained  for  a 
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homogeneous  medium  under  the  same  conditions.  Hence,  a  constant  macroscopic 
concentration  gradient  G  is  introduced  as. 

G  =  — ix.  (4.9) 

X2  —  X\ 

where  ix  is  a  unit  vector  pointing  in  the  positive  x  direction.  The  constant  macro¬ 
scopic  concentration  gradient  and  the  periodicity  of  the  medium  inevitably  lead  to 
a  fundamental  hypothesis  that  the  local  concentration  gradient  VC  is  also  spatially 
periodic  with  the  same  periodicity  as  the  porous  medium  (Adler  1992).  i.e.. 

VC(R)  =  VC(R  +  Rn),  (4.10) 


(4.15) 
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In  this  case,  jx  and  jy  are  evaluated  from  the  discrete  SPH  fluid  particles  as, 

9Ca  7Tia 
"  Ox  pa 

jx  —  x— ^  mn 


"('a 
Pa 


(4.16) 


and. 


E 


jy  =  -<k~ 


dCg  ma 
dy  pa 
ma  ’ 

Pa 


(4.17) 


where  and  are  the  local  concentration  gradients  in  x  and  y  directions. 

ox  oy 

respectively,  at  particle  a.  Values  of  dxx  and  dyx  are  obtained  from  Equations  4.14  to 
4.17.  If  another  concentration  gradient  ^  is  imposed  in  the  y  direction  only  (Figure 
4.2(b)).  the  same  procedure  leads  to  values  of  dxy  and  dyy. 

In  SPH,  the  local  concentration  gradient  VCa  is  evaluated  as, 

VCa  =  Y  —CbaVaWab.  (4.1S) 

*  pb 

If  particle  a  is  close  to  a  solid  inclusion,  the  absence  of  particle  concentrations  within 
the  inclusion  has  to  be  taken  into  account.  To  do  this.  Equation  4. IS  is  used  to  obtain 
an  intermediate  estimate  (VCQ)*, 

(VCJ*  =  Y  —CbaVaWab,  (4.19) 

t  pb 

which  is  then  corrected  by  a  factor  Ca. 

VCg  =  (VCgY/Cg,  (4.20) 


where. 


C.  =  T,  (4.21) 

Ca  reflects  the  local  number  density  of  fluid  particles  contributing  to  the  concentra¬ 
tion  at  particle  a.  Similar  approaches  have  been  used  to  improve  the  accuracy  of  first 
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derivatives  (Randles  and  Libersky  1996)  and  calculation  of  surface  tension  force  (Mor¬ 
ris  1999)  and  can  be  derived  within  the  formalism  of  element  free  Galerkin  methods 
(Belytschko  et  al.  1996). 

The  diffusion  coefficients  of  a  porous  medium  depend  on  both  medium  properties 
and  the  molecular  diffusion  coefficient  of  the  interstitial  fluid  do-  The  influence  of  do 
is  separated  out  by  defining  the  nondimensional  diffusivity  d*  (Equation  2.5S)  and 
d’  (Equation  2.61)  for  isotropic  and  anisotropic  porous  media,  respectively.  These 
values  are  functions  of  the  media  only. 

4.3  Model  Initialization  and  Execution 

4.3.1  Initialization 

For  each  simulation,  fluid  and  boundary  particles  are  initialized  in  the  same  wa\ 
as  for  the  flow  problems  described  in  Chapter  3.  Since  SPH  fluid  particles  do  not 
move  for  a  pure  diffusion  simulation,  particle  velocities  are  not  computed.  This  also 
results  in  a  constant  density  field  and  the  same  neighboring  particles  foi  evei\  SPH 
particle  throughout  a  simulation. 

After  positioning  fluid  and  boundary  particles,  linked  lists  are  created  foi  the 
particles  and  the  density  field  is  recalculated.  The  linked  lists  and  density  field  remain 
unchanged  thereafter  (i.e. ,  £2  =  x°  and  p2  =  p°  in  Equation  4.6).  For  the  pioblem  in 
Figure  4.2(a),  concentrations  on  the  left  and  right  edges  of  the  unit  cell  are  assigned 
to  be  1  and  0.  respectively,  resulting  in  a  concentration  gradient  of  For  the 

problem  in  Figure  4.2(b),  concentrations  on  the  bottom  and  top  edges  of  the  unit 
cell  are  assigned  to  be  1  and  0,  respectively,  resulting  in  a  concentration  gradient  of 
-r- .  To  expedite  the  solution,  the  concentration  field  in  the  computational  domain  is 
initialized  according  to  the  imposed  macroscopic  concentration  gradient. 
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4.3.2  Execution 

Once  the  concentration  field  is  initialized,  rates  of  change  of  particle  concentrations 
are  computed  using  Equation  4.2.  The  quintic  spline  kernel  is  employed  and  the 
smoothing  length  h  is  chosen  equal  to  2Ax.  The  time  step  is  limited  according 
to  Equation  4.4.  Concentrations  are  updated  using  the  predictor-corrector  method 
described  by  Equations  4.5  to  4.8.  The  impenetrable  solid  boundary  condition  is 
implemented  by  allowing  the  exchange  of  solute  mass  only  between  fluid  particles. 

Image  particles  and  solid  inclusions  are  created  within  a  distance  of  3h  (i.e..  one 
quintic  kernel  radius)  from  the  unit  cell  to  provide  the  necessary  "neighbors"  for 
particles  within  the  cell  (Figure  3.5).  Concentrations  of  the  image  fluid  particles  are 
updated  according  to  Equation  4.11  to  maintain  the  imposed  concentration  gradient. 

Values  of  particle  concentrations  are  recorded  periodically  during  the  course  of  a 
simulation  to  determine  the  diffusion  coefficients  of  the  porous  medium.  Changes  in 
these  values  with  time  are  used  to  determine  whether  or  not  the  diffusion  process  has 
reached  steady  state. 

4.4  Diffusion  Model  Verification 

The  SPH  diffusion  model  has  been  tested  for  many  cases.  Simulations  using  the 
method  show  excellent  agreement  with  analytical  solutions  for  diffusion  in  aqueous 
solution  and  close  agreement  with  other  solutions  for  diffusion  through  a  regular 
lattice  of  obstacles. 

4.4.1  Simulations  of  Diffusion  in  Aqueous  Solution 

Two  cases  of  diffusion  in  a  pure  aqueous  solution  with  a  molecular  diffusion  coeffi- 
2 

dent  of  d0  =  10-10  were  simulated  using  SPH.  Figure  4.3  shows  the  geometry  for 
one-dimensional  diffusion  of  a  substance  initially  confined  to  the  region  —  h  <  x  <  +h. 
The  analytical  solution  for  the  values  of  concentration  as  a  function  of  time  and  po¬ 
sition  is  (Crank  1975), 
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C=-C0 


erf 


h  -  x\ 


2y/dot  ) 


_t  (  h  +  x 

+er/U\/*f 


(4.22) 


where  erf  is  the  error  function  defined  as. 

erf  ( 2 )  =  4-  f  exp(—rf)di].  (4.23) 

7T  2  J  o 

Figure  4.4  shows  a  comparison  between  the  normalized  concentration  (^)  profiles 
obtained  using  Equation  4.22  and  SPH  for  h  =  0.0004  m  at  several  times.  The  SPH 
solution  was  obtained  with  20  particles  spanning  the  distance  h.  The  maximum  error 
of  SPH  is  about  0.5%.  confirming  the  accuracy  of  the  approach. 


Figure  4.3:  Geometry  for  one-dimensional  diffusion  in  an  aqueous  solution. 

Figure  4.5  shows  the  geometry  for  the  diffusion  of  a  substance  initially  confined  in 
a  circular  area  of  radius  R  (i.e..  diffusion  due  to  an  instantaneous  cylindrical  source). 
The  analytical  solution  for  the  values  of  concentration  as  a  function  of  time  and 
position  is  (Crank  1975), 


where  /0  is  the  modified  Bessel  function  of  the  first  kind  of  order  zero.  Figure  4.6 
shows  a  comparison  between  the  normalized  concentration  (^P)  profiles  obtained 
using  Equation  4.24  and  SPH  for  R  =  0.0005  m  at  several  times.  The  SPH  solution 
was  obtained  with  25  particles  spanning  the  distance  R.  Again,  the  results  are  in 
excellent  agreement  with  a  maximum  error  of  about  0.5%. 


[1] 


Figure  4.4:  Normalized  concentration  versus  distance  for  one-dimensional  diffusion  in 
an  aqueous  solution. 


Figure  4.5:  Geometry  for  two-dimensional  diffusion  in  an  aqueous  solution. 
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Figure  4.6:  Normalized  concentration  versus  distance  for  an  instantaneous  cylindrical 
diffusion  source. 

4.4.2  Simulations  of  Diffusion  Through  Periodic  Lattices  of  Obstacles 

A  more  challenging  test  of  the  method  involves  diffusion  through  obstacles  ar¬ 
ranged  in  a  periodic  lattice.  Steady  state  diffusion  through  a  square  lattice  of  circular 
cylinders  was  simulated  using  SPH  for  do  —  10~10  L  =  2  mm,  and  R  =  0.5  mm. 
A  solution  for  the  same  problem  was  also  obtained  using  a  FEM  program.  The  re¬ 
sults  from  both  solutions  are  compared  in  Figure  4.7.  Another  comparison  appears 
in  Figure  4.8  to  confirm  the  ability  of  the  model  to  treat  elliptical  inclusions.  The 
latter  problem  was  solved  for  do  =  10-10  L  =  2  mm.  a  =  0.8 mm,  ^  =  2,  and 
«  =  45°.  Good  agreement  is  obtained  for  the  bulk  of  the  void  space  in  both  cases. 


132 


(a) 


^  (mm) 

(b) 

Figure  4.7:  Contour  plots  of  concentration  using  (a)  FEM,  and  (b)  SPH  for  a  square 
array  of  circular  cylinders. 
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4.4.3  Simulations  of  Diffusion  Through  Composite  Media 

If  the  solid  phase  is  assigned  a  non-zero  molecular  diffusion  coefficient  (i.e.,  the 
solid  inclusion  is  penetrable  in  terms  of  tracer  diffusion),  effective  diffusion  coefficients 
for  diffusion  through  the  resulting  composite  medium  can  be  calculated.  To  do  this. 
SPH  boundary  particles  are  treated  as  fluid  particles  having  a  different  diffusion  coef¬ 
ficient.  In  this  case,  boundary  particles  must  span  the  entirety  of  solid  inclusion,  not 
only  an  annular  region  from  the  solid  surface.  The  problem  of  calculating  composite- 
media  diffusion  coefficients  dates  back  to  the  classical  studies  of  Maxwell  (1S<3)  and 
Rayleigh  (1892).  Perrins  et  al.  (1979)  furnished  a  complete  theoretical  solution  of  the 
effective  conductivity  problem  for  circular  cylinders  in  square  and  hexagonal  arrays 
by  extending  the  method  of  Rayleigh  (1892).  Table  4.1  shows  the  comparison  Ice- 
tween  the  values  of  nondimensional  diffusivity  d *  obtained  using  SPH  and  published 
by  Perrins  et  al.  (1979)  for  square  arrays  of  circular  cylinders  (Figure  4.9).  The 
solutions  are  in  close  agreement  with  a  maximum  difference  of  about  5.5CT.  The  SPH 
solutions  approach  the  values  of  Perrins  et  al.  (1979)  as  the  resolution  (i.e.,  number 
of  particles  Npart)  increases. 


Figure  4.9:  Unit  cell  of  a  square  array  of  circular  cylinders.  The  interstitial  phase  has 
a  diffusion  coefficient  of  d0  and  the  cylinder  has  a  diffusion  coefficient  of  /3do- 
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Table  4.1:  Values  of  d *  for  composite  media. 


porosity 

0 

Impart 

SPH  solution 

Perrins  et  al.  (1979) 

difference 

0.6 

2 

11600 

1.3136 

1.3080 

+0.43% 

0.6 

2 

46000 

1.3120 

1.3080 

+0.31% 

0.6 

5 

11600 

1.7785 

1.7307 

+2.76% 

0.6 

5 

46000 

1.7568 

1.7307 

+1.51  V 

0.6 

10 

11600 

2.0900 

1.9806 

+5.52% 

0.6 

10 

46000 

2.0370 

1.9806 

+2.85% 

0.4 

2 

46000 

1.5076 

1.5028 

+0.32% 

0.4 

5 

46000 

2.4251 

2.3744 

+2.14% 

0.4 

10 

46000 

3.1796 

3.0372 

+4.69% 

4.5  Diffusion  Model  Application 

The  steady  state  diffusion  problem  was  solved  using  SPH  for  spatially  periodic 
porous  media  with  circular  and  elliptical  cylinders  arranged  in  square,  staggered,  and 
hexagonal  arrays.  Diffusion  coefficients  were  obtained  from  the  steady  state  concen¬ 
tration  fields,  and  then  used  to  calculate  values  of  nondimensional  diffusivity  d*.  The 
effects  of  specific  surface  area  So,  porosity  n,  anisotropy,  and  cylinder  arrangement 

on  the  values  of  d*  are  presented.  A  total  of  154  simulations  were  completed  with  a 

2 

molecular  diffusion  coefficient  of  the  interstitial  fluid  do  =  lb  1  — .  Impenetrable 
solid  boundary  condition  applies  for  all  simulations  in  this  section  (i.e..  0  =  0).  The 
SPH  simulation  results  are  summarized  in  Tables  4.2  and  4.3. 


Table  4.2:  Summary  of  results  for  diffusion  simulations 
with  circular  cylinders. 


simulation 

So(^) 

n 

R  (mm) 

^  part 

d’ 

square  array: 

sql 

2000 

0.9 

1.0 

6938 

0.8978 

sq2 

4000 

0.9 

0.5 

693S 

0.S97S 

sq3 

20000 

0.9 

0.1 

6938 

0.8978 

sq4 

4000 

0.8 

0.5 

3384 

0.8144 

sq5 

4000 

0.7 

0.5 

6484 

0.7536 

sq6 

2000 

0.6 

1.0 

9984 

0.6968 

sq7 

4000 

0.6 

0.5 

9984 

0.6968 

sqS 

20000 

0.6 

0.1 

9984 

0.696S 

sq9 

4000 

0.5 

0.5 

9182 

0.6339 

sqlO 

2000 

0.4 

1.0 

25806 

0.564S 

sqll 

4000 

0.4 

0.5 

25806 

0.564S 

sql  2 

20000 

0.4 

0.1 

25806 

0.564S 

i  sql3 

4000 

0.3 

0.5 

71458 

0.4363 

staggered  array: 

stl 

0.9 

1.0 

12572 

0.S976 

st2 

0.9 

0.5 

12572 

0.8976 

st3 

0.9 

0.1 

12572 

0.8976 

st4 

0.8 

0.5 

8320 

0.8155 

st5 

0.7 

0.5 

11812 

0.7527 

st6 

0.6 

1.0 

11240 

0.6916 

st7 

0.6 

0.5 

11240 

0.6916 

st8 

0.6 

0.1 

11240 

0.6916 

st9 

0.5 

0.5 

20228 

0.6355 

continued  on  next  page 

Table  4.2: 

Continue. 

simulation 

So(m) 

77 

R  (mm) 

^  part 

d* 

stlO 

2000 

0.4 

1.0 

41920 

0.5633 

st  11 

4000 

0.4 

0.5 

41920 

0.5633 

st  12 

20000 

0.4 

0.1 

41920 

0.5634 

st  13 

4000 

0.3 

0.5 

115044 

0.4350 

hexagonal  array: 

hel 

2000 

0.9 

1.0 

13760 

0.8975 

he2 

4000 

0.9 

0.5 

13760 

0.89 1 5 

he3 

20000 

0.9 

0.1 

13760 

0.8976 

he4 

4000 

0.8 

0.5 

8160 

0.8156 

he5 

4000 

0.7 

0.5 

18892 

0.7579 

he6 

2000 

0.6 

1.0 

17632 

0.7023 

he? 

4000 

0.6 

0.5 

17632 

0.7023 

he8 

20000 

0.6 

0.1 

17632 

0.7024 

he9 

4000 

0.5 

0.5 

16292 

0.6543 

he  10 

4000 

0.4 

0.5 

19824 

0.6110 

hell 

2000 

0.3 

1.0 

38884 

0.5659 

he  12 

0.3 

38884 

0.5659 

hel3 

0.3 

0.1 

38884 

0.5660 

hel4 

0.2 

0.5 

95884 

0.4942 

Table  4.3:  Continue. 
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4.5.1  Effect  of  Specific  Surface  Area 


The  specific  surface  area  So  is  based  on  the  solid  volume  of  the  porous  media. 

surface  area  of  solid  inclusions 


So  = 


volume  of  solid  inclusions 


(4.25) 


In  two  dimensions,  it  becomes. 


Sn  = 


perimeter  of  solid  inclusions 


area  of  solid  inclusions 
In  this  work,  because  onlv  uniform  solid  obstacles  are  considered. 


2 

S°=R' 


(4.26) 


(4.27) 


for  c  ircular  inclusions  and, 


So  = 


4  r-V~a 
Jo 


2  cos2  x  +  b2  sin2  xdx 


nab 


(4.2S) 


for  elliptical  inclusions,  respectively. 

The  effect  of  So  on  d *  was  studied  by  calculating  values  of  d*  for  the  three  circular 
cylinder  arrays  having  So  =  2000^,  4000^,  and  20000^.  Three  porosities  of  0.4, 
O.G.  and  0.9  were  chosen  for  the  square  and  staggered  arrays,  while  three  porosities 
of  0.3.  0.6.  and  0.9  were  chosen  for  the  hexagonal  array.  The  results  are  presented  in 
rows  sql-sq3,  sq6-sq8,  sql0-sql2,  stl-st3,  st6-st8,  stl0-stl2,  hel-he3,  he6-he8,  and 
hell-hel3  of  Table  4.2.  It  is  evident  that  d*  is  not  a  function  of  So  for  any  case.  This 
is  not  surprising  since  d*  is  dimensionless  and  is  thus  not  a  function  of  any  length 
scale. 


4.5.2  Effect  of  Porosity 

The  functional  relationship  between  d‘  and  porosity  n  was  generated  using  the 
data  in  Table  4.2  and  is  shown  in  Figure  4.10.  Values  of  d*  decrease  with  decreasing  n 
for  each  array,  d*  is  essentially  independent  of  array  type  for  n  >  0.6.  For  n  <  0.6,  the 
hexagonal  array  yields  larger  values  of  d*  while  values  for  the  other  arrays  are  nearly 
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Figure  4.10:  Nondimensional  diffusivity  versus  porosity  for  square,  staggered,  and 
hexagonal  circular  cylinder  arrays. 

indistinguishable.  Due  to  the  absence  of  diffusion  within  the  solids,  the  diffusion 
process  is  controlled  by  the  distribution  of  voids  within  the  media.  The  hexagonal 
array  has  larger  diffusion  channels  with  the  same  porosity  than  square  and  staggered 
arrays,  and  this  effect  becomes  significant  when  porosity  is  low,  resulting  in  a  larger 
diffusivity  for  hexagonal  array. 

Perkins  and  Johnston  (1963)  suggested  a  single  value  of  0.7  for  d*  of  consolidated 
granular  media.  This  value  corresponds  to  a  medium  of  a  porosity  of  about  0.6 
according  to  Figure  4.10.  Fried  and  Combarnous  (19<  1)  reported  values  of  0.4  to  0.8 
for  d",  which  is  consistent  with  Figure  4.10.  Apparently,  a  value  of  g  for  d *  predicted 
bv  Saffman’s  model  (Equations  2.87  and  2.88)  only  applies  to  very  low  porosity  media. 

4.5.3  Effect  of  Anisotropy 

The  effect  of  anisotropy  was  studied  using  elliptical  inclusions  (Figures  3.11  to 
3.13).  In  Table  4.3.  the  values  of  d;x,  d*n,  d*xy ,  d*x,  as  well  as  the  principal  nondimen- 
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sional  diffusivities  d\  and  d2,  and  the  angle  9  (Figure  2.8)  are  reported.  Equations  2.49 
and  2.50  were  employed  to  obtain  d\,  d*2.  and  6.  The  values  confirm  that  d*y  =  d*yx 
in  all  cases. 

d*  d* 

Figure  4.11  shows  ^  ancj  E2.  versus  the  aspect  ratio  respectively,  for  elliptical 
cylinder  arrays  having  different  a.  d’  is  the  corresponding  nondimensional  diffusivity 
of  the  anisotropic  medium's  isotropic  counterpart,  i.e..  the  value  for  a  corresponding 
circular  cylinder  array  with  the  same  porosity  n  and  specific  surface  area  So-  Figure 
4.11  shows  that  the  anisotropy  of  the  media  increases  with  increasing  ^  and  that  the 

d*  ^ * 

values  of  ^  and  are  less  dependent  on  a  and  array  type.  Figure  4.12  presents 

d*  (j* 

and  versus  q,  respectively,  for  elliptical  cylinder  arrays  having  different  In 

d* 

Figure  4.12,  values  of  ^  and  are  grouped  in  distinguishable  bands  according  to  |. 

Figure  4.12  also  shows  that  the  effect  of  a  increases  with  increasing  It  is  concluded 

that  a  influences  the  directional  properties  of  the  anisotropy;  while  the  magnitude 

of  the  anisotropy  of  the  media,  which  can  be  defined  as  the  value  of  ^4,  is  mainly 

determined  by  the  aspect  ratio  of  the  solid  obstacles. 

Figures  4.11  and  4.12  were  generated  for  n  =  0.8.  For  the  staggered  array,  some 

d*  d* 

simulations  were  completed  for  n  —  0.5.  Figure  4.13  shows  ^  and  ^  versus 
respectively,  for  the  two  sets  of  simulations  for  staggered  array  with  porosities  of  0.5 
and  O.S.  As  can  be  seen,  a  lower  value  of  porosity  not  only  enhances  the  magnitude 
of  anisotropy  of  the  media,  it  also  increases  the  effect  of  a. 

The  relationship  between  6  and  a  is  shown  in  Figure  4.14  for  n  =  0.8.  8  is 
almost  identical  with  a  for  hexagonal  arrays,  while  is  slightly  greater  than  a  for 
staggered  arrays  and  slightly  less  than  a  for  square  arrays  when  a  is  less  than  45°. 
In  practice,  a  linear  relationship  between  9  and  a  could  be  assumed,  i.e.,  the  first 
principal  diffusion  axis  is  aligned  with  the  preferred  orientation  of  the  elongated  solid 
obstacles.  Unfortunately,  this  is  not  true  when  Figure  4.15  is  plotted,  which  also 
shows  8  versus  a  for  the  simulations  for  n  =  0.5.  A  lower  value  of  porosity  shifts  the 
first  principal  diffusion  axis  away  from  the  direction  of  particle  orientation;  however, 
the  degree  of  shifting  decreases  with  increasing  aspect  ratio 
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Figure  4.14:  6  versus  a  for  elliptical  cylinder  arrays  having  n  =  0.8. 

4.5.4  Effect  of  Cylinder  Arrangement 

The  effect  of  cylinder  arrangement  is  manifested  in  the  above  figures.  Figure 
4.10  shows  that  a  hexagonal  array  will  likely  result  in  a  higher  d*  value  than  square 
and  staggered  arrays  for  the  same  porosity.  The  effect  of  cylinder  arrangement  on 
anisotropy  is  minor  as  seen  in  Figure  4.12.  The  values  for  different  arrangements  are 
grouped  in  a  narrow  band  for  different  aspect  ratios.  It  is  interesting  to  see  that, 
hexagonal  array  always  has  a  value  of  anisotropy  in  between  those  for  the  square  and 
staggered  arrays;  and  square  array  shows  the  most  anisotropy  and  the  staggered  array 
the  least  when  a  is  less  than  about  27°. 

4.6  Summary 

Smoothed  Particle  Hydrodynamics  (SPH)  has  been  implemented  to  model  tracer 
diffusion  through  porous  media.  Comparative  studies  confirm  the  approach  is  accu¬ 
rate.  Nondimensional  diffusivities  d*  were  calculated  using  the  model  for  spatially 


Figure  4.15: 


-e —  ST  n  =  0.8  a/b  =  1.5 
■*-STn=  0.8  a/b  =  2.0 
■a —  ST  n=  0.8  a/b  =  2.5 
-e —  ST  n  =  0.5  a/b  =  1.5 
-i —  ST  n  =  0.5  a/b  =  2.0 
-* —  ST  n  =  0.5  a/b  =  2.5 


6  versus  a  for  staggered  elliptical  cylinder  arrays  having  n  —  0.5  and 
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periodic  porous  media  with  circular  and  elliptical  cylinders  arranged  in  square,  stag¬ 
gered.  and  hexagonal  arrays.  The  results  provide  guidance  for  determining  the  values 
of  d*  of  real  porous  media.  While  specific  surface  area  So  does  not  affect  d*.  porosity 
n  and  array  type  were  found  to  be  the  most  and  least  important  parameters,  respec¬ 
tively.  which  influence  the  values  of  d*.  Anisotropy  of  a  medium  is  mainly  determined 
by  aspect  ratio  ^  of  its  solid  inclusions. 
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CHAPTER  5  PORE-SCALE  TRACER  CONVECTION  AND 
HYDRODYNAMIC  DISPERSION  SIMULATIONS 

This  chapter  extends  the  SPH  numerical  models  described  in  Chapters  3  and  4 
to  study  tracer  convection  and  hydrodynamic  dispersion  in  spatially  periodic  porous 
media.  SPH  is  used  to  solve  the  Taylor  dispersion  problem  and  explore  the  charac¬ 
terization  of  dispersion  as  an  asymptotic  Fickian  process.  Discrete  SPH  particle  data 
are  analyzed  using  the  method  of  moments  and  the  advantages  and  limitations  of  the 
numerical  model  are  discussed. 

5.1  Pore-Scale  Tracer  Convection  Model 

In  the  problem  of  tracer  convection  (i.e.,  without  diffusion),  some  fluid  particles 
are  tagged  as  tracers  and  their  movement  through  void  system  of  the  porous  medium 
is  monitored.  Being  a  Lagrangian  technique,  SPH  tracks  the  trajectories  of  fluid 
particles  naturally  and  the  SPH  flow  model  (Chapter  3)  was  easily  extended  to  study 
tracer  convection.  Implementation  of  the  convection-diffusion  equation  (Chapter  4) 
was  not  needed  for  the  tracer  convection  model  because  there  is  no  mass  exchange 
between  tracer  and  tracer  or  tracer  and  carrier  fluid.  Thus,  the  solute  mass  associated 
with  each  SPH  tracer  particle  is  constant  during  the  course  of  a  simulation.  Tracer  and 
fluid  particles  are  simply  distinguished  as  “black”  and  “white”  SPH  fluid  particles, 
respectively,  in  this  model. 

Due  to  the  periodicity  of  flow,  the  tracer  convection  problem  is  solved  within  a 
single  unit  cell.  The  flow  is  first  evolved  to  steady  state;  then  a  chosen  number  of  fluid 
(white)  particles  within  the  unit  cell  are  tagged  as  tracer  (black)  particles.  If  tracers 
are  assumed  to  have  a  concentration  C,  each  tracer  particle  z  carries  a  constant  solute 
mass  Mz  of, 
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m, 

MZ  =  C — (5.1) 

Pz 

As  the  flow  continues,  tracer  and  fluid  particles  are  wrapped  around  the  unit  cell 
in  the  same  manner  as  described  in  Chapter  3  (Figure  3.6).  However,  information 
regarding  the  wrapping  of  tracer  particles  around  the  unit  cell  is  recorded  and  used  to 
construct  the  real  tracer  field  in  the  corresponding  spatially  periodic  porous  medium. 
For  example,  in  Figure  5.1(a).  if  the  tracer  particle  is  known  to  have  crossed  the  upper 
unit  cell  boundary  once  and  the  right  unit  cell  boundary  twice,  its  real  position  in  the 
spatially  periodic  porous  medium  is  shown  in  Figure  5.1(b).  Calculations  are  needed 
for  only  one  unit  cell  to  simulate  tracer  convection  using  this  method. 

5.2  Pore-Scale  Hydrodynamic  Dispersion  Model 

With  implementations  of  both  flow  (Chapter  3)  and  convection-diffusion  (Chapter 
4)  equations,  a  pore-scale  tracer  hydrodynamic  dispersion  model  was  developed  using 
SPH.  The  SPH  model  easily  simulates  the  process  of  tracer  hydrodynamic  dispersion 
in  principle,  however,  in  practice,  simulations  are  limited  by  computation  time. 

In  spatially  periodic  porous  media,  while  the  flow  field  is  also  periodic  in  nature, 
this  is  not  true  for  the  concentration  field.  As  a  result,  in  the  problem  of  simulating 
tracer  hydrodynamic  dispersion,  a  computing  domain  consisting  of  multiple  unit  cells 
is  needed  for  evolving  the  concentration  field.  For  example,  although  the  initial  tracer 
plume  is  confined  within  one  unit  cell  in  Figure  5.2(a),  the  tracer  plume  at  a  later 
time  spans  to  six  unit  cells  in  Figure  5.2(b).  An  efficient  algorithm  to  model  the 
evolution  of  concentration  which  takes  advantage  of  the  periodic  nature  of  flow  field 
is  described  below. 

Figure  5.3  shows  the  initial  concentration  field  for  a  problem  of  tracer  dispersion. 
The  concentration  computing  domain  consists  of  six  unit  cells  which  are  identified 
by  the  (i,  j)  coordinates.  For  this  tracer  hydrodynamic  dispersion  model,  the  flow 
field  is  calculated  in  a  single  unit  cell  corresponding  to  i  —  1  and  j  —  1.  Every  SPH 


Figure  5.2:  Tracer  plumes  for  the  problem  of  tracer  hydrodynamic  dispersion:  (a) 
initial  tracer  plume,  and  (b)  tracer  plume  at  time  t. 
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fluid  particle  only  has  one  set  of  flow  related  quantities  (e.g.,  velocity  and  density) 
and  carries  as  many  concentration  values  as  the  number  of  unit  cells  which  constitute 
the  concentration  computing  domain.  In  this  case,  every  fluid  particle  has  six  values 
of  concentration,  one  for  each  unit  cell  in  the  concentration  computing  domain.  The 
concentration  field  is  determined  by  the  particle  concentration  values  for  each  unit 
cell.  However,  particle  positions  in  unit  cell  (i,  j )  other  than  unit  cell  (1.  1)  are  not 
needed  during  the  course  of  a  simulation.  According  to  Equation  4.2.  the  evolution 
of  concentration  field  is  driven  by  local  concentration  gradients  and  only  the  relative 
positions  between  particles  are  needed  to  calculate  the  rates  of  change  for  concentra¬ 
tion.  The  relative  positions  between  particles  are  calculated  from  the  position  values 
in  unit  cell  (1,  1)  and  they  are  the  same  for  every  unit  cell.  For  the  purpose  of  analyz¬ 
ing  and  visualizing  the  concentration  field,  SPH  particles  (both  fluid  and  boundary 
particles)  are  mapped  to  other  unit  cells  from  unit  cell  (1,  1)  in  the  concentration 
computing  domain. 

To  perform  the  simulation  in  Figure  5.3,  the  flow  field  is  evolved  to  steady  state 
with  all  particle  concentrations  equal  to  zero.  In  order  to  initialize  the  dispersion 
process,  particles  in  area  A  (Figure  5.4)  are  assigned  Co  =  1  for  unit  cell  (1,  1). 
Likewise,  particles  in  area  B,  C,  D  (Figure  5.4)  are  assigned  Co  —  1  for  unit  cells  (1, 
2).  (2.  2),  and  (2,  1).  respectively. 

Once  the  concentration  field  is  initialized,  the  flow  continues  according  to  the 
methods  described  in  Chapter  3  and  the  concentration  field  is  evolved  according  to 
the  methods  described  in  Chapter  4.  As  fluid  particles  are  wrapped  around  the 
unit  cell  (1,  1)  and  recalling  that  all  SPH  particles  physically  reside  in  unit  cell  (1, 
1).  concentration  values  of  the  wrapped  particles  are  reassigned  accordingly.  For 
example,  as  fluid  particle  a  (Figure  5.5)  is  wrapped  back  into  unit  cell  (1,  1),  its  value 
of  concentration  for  unit  cell  ( i ,  j )  after  wrapping  equals  to  the  value  of  concentration 
for  unit  cell  (i  -  1,  j)  before  wrapping. 

Values  of  concentration  must  also  be  correctly  assigned  for  image  fluid  particles. 
For  example,  if  particle  c  is  the  image  particle  of  particle  b  in  unit  cell  (1,  1)  (Figure 


i=  1 


i=  2 


3 


Figure  5.3:  Initial  concentration  field  for  a  problem  of  tracer  hydrodynamic  dispersion 


in  a  spatially  periodic  porous  medium. 


y 


i  =  1 


Figure  5.4:  Unit  cell  (1,1)  for  the  problem  in  Figure  5.3. 
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5.5),  the  value  of  concentration  for  particle  c  corresponding  to  unit  cell  (i,  j)  is  equal 
to  the  value  of  concentration  for  particle  b  corresponding  to  unit  cell  (i  +  l.  j). 


Figure  5.5:  Wrapping  fluid  particles  and  creating  image  particles  to  simulate  tracer 
dispersion  in  a  spatially  periodic  porous  medium. 

The  concentration  computing  domain  used  in  this  work  consists  of  a  constant 
number  of  unit  cells.  A  large  concentration  computing  domain  is  unnecessary  during 
the  initial  stage  of  a  simulation.  For  example,  the  unit  cells  corresponding  to  column 
i  =  3  in  Figure  5.3  are  not  needed  when  the  dispersion  simulation  begins.  It  is  possible 
to  dynamically  change  the  size  of  concentration  computing  domain  to  accommodate 
the  evolution  of  concentration.  However,  this  requires  dynamically  determining  the 
tracer  plume  front.  The  advantages  of  a  dynamically  changing  concentration  com¬ 
puting  domain  during  the  course  of  a  simulation  are  not  explored  in  this  work.  For  a 
constant  concentration  computing  domain,  a  simulation  terminates  when  the  plume 
front  crosses  the  boundary  of  the  domain.  A  meaningful  simulation  of  hydrodynamic 
dispersion  usually  requires  several  tens  of  unit  cells  for  modeling  the  evolution  of 
concentration.  This  significantly  lengthens  the  computation  time  for  each  simulation. 

Since  SPH  uses  an  explicit  time  integration  scheme,  the  time  step  is  limited  by 
stability  constraints.  For  a  low  Reynolds  number  flow  simulation,  the  time  step  is 
usually  limited  by  fluid  viscosity  according  to  Equation  3.47, 
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t  <  0.125  —  . 

v 

The  time  step  for  properly  integrating  particle-to-particle  diffusion  is  (Equation  4.4). 

Ir 

r  <  0.125—. 
do 

As  liquids  have  a  Schmidt  number  Sc  =  on  the  order  of  103  to  10*.  the  maximum 
time  step  is  limited  by  Equation  3.47  for  a  simulation  of  liquid-typo  tracer  dispersion. 
As  a  result,  the  chosen  time  step  is  relatively  small  for  evolving  concentrations  and 
lengthy  computer  simulations  are  needed  to  simulate  dispersion  phenomena,  due  to 
the  large  concentration  computing  domain. 

5.3  Method  of  Moments 


Spatial  moments  of  tracer  distribution  provide  an  estimation  of  the  permeability 
and  dispersivity  of  porous  media  (Aris  1956;  Horn  1971;  Brenner  1980a).  In  the 
SPH  tracer  convection  model,  the  tracer  field  is  characterized  by  individual  tracer 
particle  position  R  and  its  associated  solute  mass  Mz.  The  zeroth  moment  of  tracer 
distribution  is. 


dM 
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J  dMo  MZfi 


(5.2) 


Mo  will  have  a  constant  value  of  unity  since  AT 
moment  of  tracer  distribution  is. 


AT,o  in  this  model.  The  first 
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Mi  represents  the  position  of  the  center  of  tracer  mass  and  is  used  to  calculate  the 
seepage  velocity  \s  as, 
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The  centered  second  moment,  is, 


M,  = 


J (R  -  Mi)2dM  E(R-' 
[  dM0 


represents  the  spreading  of  a  tracer  field  about  its  center  and  is  used  to  calculate 
the  dispersion  coefficient  as, 


Equation  5.6  is  consistent  with  Equation  2.150. 

For  the  SPH  tracer  hydrodynamic  dispersion  model,  the  solute  mass  of  a  given 
particle  is  not  constant  and  Equations  5.2.  5.3,  and  5.5  become. 
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where  summations  are  performed  over  all  fluid  particles.  Due  to  molecular  diffusion, 
every  SPH  fluid  particle  could  potentially  carry  some  solute  mass  in  the  tracer  hydro- 
dynamic  dispersion  model.  The  zeroth  moment  (Equation  5.7)  will  have  a  constant 
value  of  unity  if  tracer  mass  is  conserved  during  a  simulation  and  this  provides  a  first 
check  of  the  model. 

If  tracer  dispersion  is  a  Fickian  process,  the  dispersive  flux  is  proportional  to  the 
concentration  gradient  (Equation  2.67).  It  is  usually  believed  that  tracer  dispersion 
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is  an  asymptotic  Fickian  process.  Under  such  conditions,  the  dispersive  flux  can  be 
characterized  by  a  constant  dispersion  coefficient  after  a  characteristic  time  tc.  which 
is  the  time  required  for  the  tracer  to  sample  interstitial  space  of  the  unit  cells  (Brenner 
1980a).  However,  there  are  some  studies  in  literature  where  dispersivity  may  never 
approach  an  asymptotic  value  (section  2.3.7). 

5.4  Interpretation  of  Tracer  Diffusion  Using  Method  of  Moments 

Figure  5.6  shows  the  geometry  for  one-dimensional  diffusion  of  a  substance  initially 

confined  to  a  finite  region.  The  solid  configuration  is  a  periodic  square  array  of 

2 

circular  cylinders.  The  problem  was  solved  for  do  —  lO-10^-.  R  —  0.5  mm.  n  = 
O.S.  and  h  =  0.6  mm.  The  method  of  moments  was  used  to  calculate  the  diffusion 


Solute  is  confined  to  this  region  at  t  =  0 

Figure  5.6:  Initial  condition  for  a  tracer  diffusion  through  a  square  array  of  circular 
cylinders. 

coefficient  in  the  x  direction  dx  using  Equation  5.6  with  A t  =  1570  s.  dx  was  then 
used  to  calculate  the  nondimensional  diffusivity  in  the  x  direction  d*x  (Equation  2.58). 
Figure  5.7  shows  d*x  as  a  function  of  t.  Values  of  d*x  approach  an  asymptotic  value 
of  0.S13S  with  time.  For  the  same  media  configuration,  dx  was  found  to  be  0.8144 
using  the  steady  state  diffusion  approach  presented  in  Chapter  4  (Table  4.2).  The 
method  of  moments  produced  a  result  which  is  consistent  with  that  obtained  using  the 
steady  state  diffusion  approach.  However,  the  method  of  moments  simulation  used  a 
concentration  computing  domain  consisting  of  30  unit  cells  and  therefore  required  a 
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f  (x  103  s) 


Figure  5.7:  Nondimensional  diffusivity  versus  time  for  diffusion  through  a  square 
array  of  circular  cylinders. 

much  longer  time  to  run. 

5.5  Simulations  of  Taylor  Dispersion 

Taylor  dispersion  between  stationary  infinite  plates  at  y  =  0  and  y  =  L  (Figure 
3.9)  was  simulated  using  the  tracer  hydrodynamic  dispersion  model.  The  method  of 
moments  was  employed  to  yield  Taylor  dispersion  coefficients  which  were  compared 
with  values  from  Taylor's  analytical  solution. 

Taylor  (1953)  derived  Equation  2.85  to  describe  the  average  cross  section  concen¬ 
tration  in  a  cylindrical  straight  tube.  Equation  2.85  indicates  that  such  dispersion 
may  be  characterized  as  an  asymptotic  diffusional  process  with  an  effective  disper- 

r>2  2 

sion  coefficient  DTayior  =  d0  +  For  the  system  shown  in  Figure  3.9;  the  Taylor 

dispersion  coefficient  becomes  (Horne  and  Rodriguez  1983;  Hull  et  al.  1987), 


D 


Taylor 


Z,V 

=  d0  + 


210  <V 


(5.10) 


162 


where  vs  is  the  average  channel  velocity  determined  by, 

Vs  =  t \~FL2. 

12  v 

The  Reynolds  number  and  Peclet  number  are. 

LV' 


Re 


and. 


Pe  = 


L .v, 
do  ’ 


(5.11) 


(5.12) 


(5.13) 


respectively.  The  relationship  between  Drayior  and  Pe  is, 

'Taylor  _  i  i _ —Pe2 

d0  210 


(5.14) 


SPH  was  used  to  simulate  Taylor  dispersion  between  stationary  infinite  plates 
with  L  =  0.001m.  p  =  103-H.  v  =  1CT6  TJf-.  d0  =  lO^6  F  =  0.05.  0.10.  and 

77T  5  S 

0.15  Rj.  and  50  fluid  particles  spanning  the  distance  L.  Poise.uille  flow  was  evolved  in 
a  square  unit  cell  of  the  size  L  x  L  and  the  size  of  the  dispersion  computing  domain 
was  40L  x  L  (Figure  5.8). 


Figure  5.8:  Initial  condition  ( t  =  0)  for  Taylor  dispersion  between  two  infinite  plates. 

2 

A  Schmidt  number  of  1  (i.e.,  v  =  do  =  10-6  ^-)  was  assumed  for  the  fluid  so 
that  the  time  step  due  to  viscosity  would  be  comparable  with  that  due  to  molecular 
diffusion.  While  the  chosen  value  of  u  is  reasonable  for  a  liquid,  the  value  of  do  is 
more  appropriate  for  a  gas.  For  the  simulations  considered  here,  the  maximum  time 


163 


step  is  At  =  0.0002  s  as  determined  by  Equation  3.47.  To  observe  the  asymptotic 
Taylor  dispersion  behavior,  time  t  must  evolve  to  a  value  greater  than  tc ,  which  is 
defined  appropriately  for  this  problem  geometry  as. 


according  to  Equation  2.84.  Equation  5.15  yields  tc  =  0.25  s  for  d{)  =  10-6  A 
simulation  with  a  value  of  d0  appropriate  for  a  liquid,  such  as  d0  =  10“ 10  could  be 
conducted:  however,  the  value  of  tc  would  then  be  2500  s  and  at  least  12.500.000  time 
steps  (At  =  0.0002  s)  would  be  required  for  t.  >  tc.  Furthermore,  as  the  tracer  plume 
would  be  transported  much  further  with  increasing  time,  a  much  larger  dispersion 
computing  domain  would  be  needed  to  contain  the  tracer  body.  Within  practical 
limits,  such  a  simulation  could  not  be  completed  using  currently  available  computers. 

After  the  Poiseuille  flow  was  evolved  to  steady  state,  the  Taylor  dispersion  simu¬ 
lation  started  with  the  initial  condition  shown  in  Figure  5.8.  Figure  5.9  shows  plots 
of  the  tracer  concentration  field  at  four  times  for  a  simulation  with  F  =  0.10  The 
plots  were  generated  directly  using  discrete  SPH  data  with  the  grey-scale  representing 
the  value  of  concentration. 

Tracer  dispersion  was  analyzed  using  the  method  of  moments.  Figure  5.10  is  a 

plot  of  the  zeroth  moment  Mo  as  a  function  of  time.  The  plot  shows  that  tracer  mass 

was  conserved  in  the  simulations.  Figure  5.11  is  a  plot  of  the  first  moment  in  the  x 

direction  M\x  as  a  function  of  time.  Values  of  M\x  were  computed  by  subtracting 

the  first  moment  at  time  t  =  0  from  that  at  time  t.  Figure  5.11  shows  that  the 

center  of  tracer  mass  moved  at  a  constant  velocity  in  each  case.  The  average  channel 

velocity  vs  was  Calculated  using  Equation  5.4  with  At  =  0.023,  0.019  and  0.020  s  for 

F  =  0.05,  0.10,  and  0.15™,  respectively,  and  is  shown  in  Figure  5.12.  Figure  5.12 

s 

indicates  that  the  flow  fields  were  numerically  stable  after  steady  state  conditions 
were  reached.  SPH  solutions  for  vs  differ  from  the  analytical  solutions  by  a  maximum 


of  1.2%. 
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Figure  5.13  is  a  plot  of  centered  second  moment  in  the  x  direction  M2x  as  a  function 

of  time.  Values  of  M2x  were  computed  by  subtracting  centered  second  moment  at 

time  t  =  tc  from  that  at  time  t.  Figure  5.13  shows  that,  in  each  case,  tracer  spreading 

about  its  center  grows  linearly  with  time  for  t  >  tc.  The  dispersion  coefficient  in 

the  x  direction  Dx  was  calculated  using  Equation  5.6  with  At  =  0.023.  0.019  and 

0.020 .s  for  F  =  0.05.  0.10.  and  0.15  respectively,  and  is  shown  in  Figure  5.14.  It 

s 

is  seen  that  for  t  >  tc  =  0.25s.  Dx  fluctuates  about  an  average  value.  Asymptotic 
values  of  Dx  were  taken  as  Taylor  dispersion  coefficients  Djayior  and  were  found  to  be 
1.0605  x  10"6.  1.3034  x  10~6.  and  1.684S  x  10"6  ^  for  F  =  0.05.  0.10.  and  0.15  W,. 
respectively.  SPH  solutions  for  DTayior  were  less  than  the  analytical  solutions  by  a 
maximum  of  3.4%.  Taylor  dispersion  is  essentially  a  one-dimensional  phenomenon. 
The  changes  in  the  first  moment  and  centered  second  moment  in  the  y  direction  were 
found  to  have  a  maximum  magnitude  of  10~8m  and  10-9m2,  respectively,  for  SPH 
simulations.  These  values  likely  result  from  numerical  errors.  A  summary  of  results 
for  Taylor  dispersion  simulations  is  shown  in  Table  5.1. 
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1 .0002 


Figure  5.10:  The  zeroth  moment  of  tracer  distribution  versus  time  for  Taylor  disper¬ 
sion  between  two  infinite  plates. 
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Figure  5.11:  The  first  moment  in  the  x  direction  of  tracer  distribution  versus  time 
for  Taylor  dispersion  between  two  infinite  plates. 
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Figure  5.12:  Average  channel  velocity  versus  time  for  Taylor  dispersion  between  two 
infinite  plates. 
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Figure  5.13:  Centered  second  moment  in  the  x  din 
time  for  Taylor  dispersion  between  two  infinite  pla 
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Figure  5.14:  Hydrodynamic  dispersion  coefficient  in  the  x  direction  versus  time  for 
Taylor  dispersion  between  two  infinite  plates. 
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5.6  Simulations  of  Tracer  Convection 

In  this  section,  the  process  of  tracer  convection  (i.e..  without  diffusion)  in  two- 
dimensional  spatially  periodic  porous  media  is  studied.  Five  simulations  were  com¬ 
pleted  for  tracer  convection  through  spatially  periodic  porous  media  with  circular 
cylinders  arranged  in  square,  staggered,  and  hexagonal  arrays.  The  fluid  was  mod¬ 
eled  as  water,  i.e..  p  =  103  H  and  //  =  1<T6  ^  and  tracer  convection  was  driven 

m 

by  body  force  F.  As  shown  in  Figure  5.15.  7  is  the  angle  between  the  direction  of  F 
and  the  positive  x  direction.  The  direction  of  F  is  denoted  by  L.  while  the  direction 
perpendicular  to  L  is  denoted  by  T.  Table  5.2  presents  a  summary  of  information 
for  the  tracer  convection  simulations.  Figures  5.16  to  5.20  show  the  tracer  fields  at 
four  times  for  each  simulation.  The  tracer  fields  were  analyzed  using  the  method  of 
moments  for  seepage  velocity  vs,  effective  porosity  ntj gr,  and  mechanical  dispersion  co¬ 
efficients  DmL  and  DmT ■  Tortuosity  T  of  the  media  was  also  calculated.  A  summary 
of  results  for  the  tracer  convection  simulation  analysis  is  presented  in  Table  5.3. 


Figure  5.15:  Body  force  F  direction  L  and  perpendicular  direction  T. 


Table  5.2:  Summary  of  information  for  tracer  convection 
simulations. 


Table  5.3:  Summary  of  results  for  tracer  convection  sim¬ 
ulation  analysis. 


1  1.051  9.30  x  10"5  0.047  7.35  x  10~5  0.790  0.8 

2  1.139  9.06  x  10"5  0.045  7.17  x  10"5  0.791  0.S 

3  1.067  9.08  x  10-5  0.045  7.18  x  10"5  0.791  0.8 

4  1.094  6.57  x  10~5  0.033  5.19  x  10~5  0.790  0.8 

5  1.082  6.56  x  10~5  0.033  5.18  x  10“5  0.790  0.S 
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Figure  5.18:  Tracer  convection  through  a  staggered  array  of  circular  cylinders  (7 
30°).  Tracer  particles  are  shown  in  black. 
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Figure  5.20:  Continue. 
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5.6.1  Tortuosity  of  Porous  Media 

In  SPH,  trajectories  of  individual  fluid  masses  (SPH  tracer  particles)  can  be 
recorded  as  they  travel  through  the  void  system  of  the  medium.  Figures  5.22  to 
5.26  show  the  trajectories  of  four  tracer  particles  for  each  simulation  in  Table  5.2. 
Tortuosity  of  the  porous  medium  T  was  calculated  as. 

Lg  2 

~Z7 

T=  — - ,  (3.16) 

AT 

+  ' tracer 

where  Le%z  and  Lz  are  the  length  of  trajectory  and  the  straight-line  distance  traveled 
by  tracer  particle  2,  respectively  (Figure  5.21).  and  Ntracer  is  the  total  number  of 
tracer  particles  in  the  flow  field. 


Uy 

Unit  cell  z 

Figure  5.21:  Trajectory  Lej2  and  straight-line  distance  L  traveled  by  tracer  particle 


Figure  5.24:  Trajectories  of 
(l  =  30°). 
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Figure  5.26:  Trajectories  of  tracer  particles  in  a  hexagonal  array  of  circular  cylinders 
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Figure  5.27  shows  the  tortuosity  of  the  media  as  a  function  of  time  for  each 
tracer  convection  simulation  in  Table  5.2.  A  constant  tortuosity  value  exists  for 
each  case,  which  is  presented  in  Table  5.3.  A  few  interesting  observations  can  be 
made  from  Figure  5.27.  Tortuosity  values  are  considerably  less  than  the  commonly- 


1.154 


Figure  5.27:  Tortuosity  versus  time  for  simulations  in  Table  5.2. 

used  value  of  \/2  =  1.414,  which  was  originally  proposed  by  Carman  (1937)  for 
unconsolidated  porous  aggregates.  It  is  also  apparent  that  tortuosity  is  a  function 
of  media  geometry.  For  the  three  simulations  with  7  =  0,  the  square  and  staggered 
arrays  have  the  smallest  and  largest  tortuosity,  respectively.  For  the  square  array 
with  7  =  0,  most  of  the  tracer  particles  will  have  an  unobstructed  travel  path  (Figure 
5.22),  which  produces  a  tortuosity  value  close  to  unity.  The  staggered  array  introduces 
more  obstacles  into  the  flow  field  and  tortuosity  increases  as  a  result.  Although  the 
hexagonal  array  has  the  same  solid  surface  area  in  the  unit  cell  as  the  staggered  array, 
it  has  larger  flow  channels  which  results  in  a  tortuosity  larger  than  the  square  array, 
but  less  than  the  staggered  array. 

Not  only  a  function  of  media  geometry,  tortuosity  is  also  a  function  of  the  direction 
of  the  applied  body  force.  The  effect  of  body  force  direction  on  tortuosity  is  apparent 
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for  the  staggered  array  (simulations  2  and  3  in  Table  5.2).  For  7  -  30°,  fewer  tracer 
particles  must  travel  around  the  center  solid  in  the  staggered  array  than  that  if  7  =  0. 
This  results  in  a  lower  tortuosity  value  for  7  =  30°.  Similar  observations  can  be  made 
regarding  the  effect  of  7  for  the  hexagonal  array. 

5.6.2  Seepage  Velocity  and  Effective  Porosity  of  Porous  Media 

Seepage  velocity  for  each  tracer  convection  simulation  was  calculated  using  Equa¬ 
tion  5.4.  Figure  5. 28  shows  seepage  velocity  in  the  L  direction  vsl  and  in  the  T 
direction  vsj-  as  a  function  of  time  for  simulation  1  in  Table  5.2.  vsi  fluctuates  about 
an  average  value  of  vs  =  9.30  x  10"5  ™  within  a  5%  range,  which  indicates  that 
the  tracer  flow  field  is  numerically  stable.  vst  is  essentially  zero  as  expected.  The 
fluctuations  of  vsT  about  zero  likely  result  from  numerical  errors.  Other  values  of  vsL 
and  v$t  for  other  simulations  in  Table  5.2  show  similar  patterns  when  plotted  as  a 
function  of  time.  Values  of  vs  are  presented  in  Table  5.3.  which  were  also  used  to 
calculate  Re  =  for  the  simulations.  The  small  values  of  Re  indicate  creeping 
flow  conditions. 

Effective  porosity  neg  of  the  medium  was  calculated  as. 

neff  =  (5-17) 

Vs 

where  v  is  the  Darcy  velocity  in  body  force  direction  calculated  using  Equations  3.62 
and  3.63  for  the  steady  state  flow  field.  Values  of  neg  are  presented  in  Table  5.3. 
The  effective  porositv  excludes  “dead  end  ’  pores  and  thus  represents  the  fraction  of 
pore  space  that  is  available  to  transmit  fluids.  For  the  media  type  considered  in  this 
work.  neff  should  be  equal  to  n.  The  small  difference  between  neS  and  n  in  Table  5.3 
is  due  to  the  placement  of  SPH  boundary  particles.  As  the  first  layer  of  boundary 
particles  is  placed  on  the  solid  surface,  the  porous  medium  simulated  actually  has  a 
porosity  slightly  smaller  than  the  nominal  value.  Although  neg  is  nearly  equal  to  n 
for  uniform  granular  materials,  it  is  significantly  less  than  n  for  clays  due  to  their 
more  complex  microfabric  of  clay  particle  clusters  and  intercluster  voids. 
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Figure  5.28:  Seepage  velocity  versus  time  for  simulation  1  in  Table  5.2. 

5.6.3  Mechanical  Dispersion  Coefficient  of  Porous  Media 

Mechanical  dispersion  coefficients  for  tracer  convection  were  calculated  using  Equa¬ 
tion  5.6.  Figures  5.29  to  5.33  show  -ATr-  Dmi,  and  Dmr  as  a  function  of  time 
for  each  tracer  convection  simulation.  A/2L-  A/2r,  Dmi ,  and  DmT  denote  centered 
second  moments  and  mechanical  dispersion  coefficients  in  the  L  and  T  directions, 
respectively.  Values  of  M2i  and  M2t  were  computed  by  subtracting  centered  second 
moments  at  time  t  =  0  from  those  at  time  t. 

M2l  and  Dmi  show  different  types  of  behavior  depending  on  the  body  force  di¬ 
rection.  When  F  is  aligned  with  a  line  of  symmetry  for  the  solid  inclusions,  periodic 
streamlines  exist  and  the  increase  of  M2i  with  t  is  nearly  quadratic.  As  a  result, 
Dmi  increases  nearly  linearly  with  t  (Figures  5.29,  5.30,  and  5.32).  However,  when  F 
is  not  aligned  with  a  line  of  symmetry  for  the  solid  inclusions,  irregular  streamlines 
exist  and  M2L  shows  a  more  complicated  behavior  with  t  (Figures  5.31  and  5.33). 
Although  M2L  generally  increases  with  t  in  this  case,  it  may  decrease  with  t  for  short 
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(a) 


(b) 


Figure  5.30:  (a)  Centered  second  moment,  and  (b)  mechanical  dispersion  coefficient 
versus  time  for  tracer  convection  through  a  staggered  array  of  circular  cylinders  (7  = 
0°). 
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Figure  5.32:  (a)  Centered  second  moment,  and  (b)  mechanical  dispersion  coefficient 
versus  time  for  tracer  convection  through  a  hexagonal  array  of  circular  cylinders 

(7  =  0°). 
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periods,  which  results  in  a  negative  value  of  DmL.  Figures  5.31(b)  and  5.33(b)  show 
that  DmL  does  not  grow  with  time  and  instead  fluctuates  about  an  average  value. 
When  F  is  aligned  with  a  line  of  media  symmetry.  M2T  and  Dmr  are  practically  zero 
(Figures  5.29.  5.30.  and  5.32).  When  F  is  not  aligned  with  a  line  of  media  symmetry. 
M2t  shows  small  fluctuations  about  zero  which  are  amplified  when  Dmr  is  calculated 
(Figures  5.31  and  5.33).  However,  the  magnitude  of  DmT  is  consistently  smaller  than 
that  of  DmL- 

The  relationship  between  centered  second  moment  M2  and  time  can  also  be  studied 
using  a  log-log  scale.  It  is  usually  believed  that  log{M2)  and  log(t)  have  a  linear 
relationship  as, 

log  {Mi)  =  slog(t ),  (5.18) 

where  s  provides  an  indication  of  dispersion  behavior  (Cushman  1997).  A  value  of 
,s  =  1  corresponds  to  Fickian  behavior.  Figures  5.34  to  5.3S  show  log-log  plots  for 
each  tracer  convection  simulation  with  a  line  fitted  using  least  squares.  The  slope 
of  the  fitted  line  s  and  the  regression  coefficient  R  are  denoted  on  the  plots.  An 
approximately  linear  relationship  exists  for  log(M2i)  with  log(t)  (R  >  0.92).  however, 
it  is  not  appropriate  to  fit  a  line  for  the  variation  of  log(M2r)  versus  log(t )  (R  <  0.1). 
The  linear  fits  for  M2l  are  quite  good  if  F  is  aligned  with  a  line  of  media  symmetry 
(Figures  5.34  and  5.37).  log-log  plots  cannot  be  made  for  M2t  since  these  values  are 
practically  zero.  When  F  is  not  aligned  with  a  line  of  media  symmetry,  s  has  a  value 
close  to  unity  for  the  log(M2L)-log(t)  relationship,  which  suggests  an  approximate 
Fickian  dispersion  behavior  for  this  case. 

From  the  tracer  convection  simulations  under  laminar  flow  conditions,  an  asymp¬ 
totic  mechanical  dispersion  coefficient  was  not  found  for  two-dimensional  spatially 
periodic  porous  media  when  F  is  aligned  with  a  line  of  media  symmetry,  even  for 
large  times.  It  is  necessary  to  reexamine  the  fundamentals  of  the  theory  of  Fickian 
approximation  to  understand  why  the  phenomenon  does  not  exist  for  this  case.  The 
Fickian  approximation  assumes  that  the  travel  time  for  an  individual  tracer  particle 


Figure  5.34:  log{M1L)  versus  log(t)  for  tracer  convection  through  a  square  array  of 
circular  cylinders  (7  =  0°). 


Figure  5.35:  /o^(M2t)  versus  log(t)  for  tracer  convection  through  a  staggered  array 
of  circular  cylinders  (7  =  0°). 
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Figure  5.36:  log( A/2)  versus  log{t)  for  tracer  convection  through  a  staggered  array  of 
circular  cylinders  (7  =  30°). 


Figure  5.37:  /og(M2L)  versus  log(t)  for  tracer  convection  through  a  hexagonal  array 
of  circular  cylinders  (7  =  0°). 
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log(t) 

Figure  5.38:  log(M2)  versus  log(t)  for  tracer  convection  through  a  hexagonal  array  of 
circular  cylinders  (7  =  45°). 

is  much  larger  than  the  time  interval  during  which  its  successive  velocities  are  corre¬ 
lated.  Total  displacement  is  considered  as  the  sum  of  a  large  number  of  elementary 
displacements  which  are  statistically  independent  of  one  another.  Each  particle  in 
the  tracer  cloud  is  assumed  to  transverse  all  possible  variations  in  the  velocity  field. 
These  assumptions  lead  to  a  normal  distribution  of  a  cloud  of  initially  close  particles 
according  to  the  central  limit  theorem  for  an  ergodic  system,  which  is  characterized  by 
a  constant  dispersivity  (Bear  1972).  In  a  random  porous  medium,  these  assumptions 
are  reasonable  to  some  extent. 

Constant  dispersivity  implies  that  a  single  velocity  path  is  statistically  representa¬ 
tive  of  all  velocity  paths.  This  means  that  a  tracer  particle  released  anywhere  in  the 
unit  cell  of  a  spatially  periodic  porous  medium  eventually  samples  the  entire  unit  cell 
by  convection  alone.  When  F  is  aligned  with  a  line  of  media  symmetry,  Figures  5.22, 
5.23.  and  5.25  show  that  different  subsets  of  tracer  particles  may  experience  entirely 
different  velocity  paths  in  periodic  porous  media.  Preferred  flow  paths  and  spatially 
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periodic  flow  patterns  are  apparent  in  those  figures.  In  this  case,  the  structure  of  the 
periodic  porous  medium  and,  therefore,  the  velocity  experienced  by  a  tracer  particle 
traversing  a  streamline,  remain  correlated  throughout  all  space.  The  distance  between 
two  tracer  particles  released  on  different  streamlines  increases  with  time.  As  such, 
there  is  no  mechanism  for  tracer  particles  traveling  through  periodic  porous  media 
under  this  condition  to  mix  across  streamlines  in  pure  tracer  convection  simulations. 

If  F  is  not  aligned  with  a  line  of  media  symmetry  (Figures  5.24  and  5.26).  although 
streamlines  still  cannot  cross  each  other  under  laminar  flow  conditions,  irregular  pat¬ 
terns  of  streamlines  seem  to  provide  the  possibility  for  a  tracer  particle  to  sample  the 
whole  unit,  cell  during  the  course  of  a  simulation.  This  results  in  a  mixing  process 
which  can  be  approximated  as  Fiekian.  Another  mechanism  which  causes  tracer  mass 
to  mix  across  streamlines  is  molecular  diffusion. 

5.7  Simulations  of  Tracer  Hydrodynamic  Dispersion 

Due  to  computational  limitations,  extensive  simulations  of  tracer  hydrodynamic 
dispersion  through  spatially  periodic  porous  media  are  not  presented.  As  an  asymp¬ 
totic  Fiekian  behavior  of  pure  tracer  convection  was  not  found  for  two-dimensional 
spatially  periodic  porous  media  when  F  is  aligned  with  a  line  of  media  symmetry,  the 
issue  of  whether  tracer  hydrodynamic  dispersion  under  this  condition  can  be  described 
as  an  asymptotic  Fiekian  process  is  explored  in  this  section. 

One-dimensional  tracer  hydrodynamic  dispersion  through  a  periodic  square  array 

of  circular  cylinders  (n  =  0.8  and  R  =  0.5  mm)  was  simulated  for  p  =  lO3-^-, 

v  '  m 

v  =  10"6  F  =  0.001  and  d0  =  10~10  10~9  and  10~8  The  body 

s 

force  F  was  applied  in  the  positive  x  direction,  i.e.,  7  =  0°,  and  the  induced  steady 
state  seepage  velocity  vs  was  9.30  x  10-5  FI,  For  the  simulations.  Re  =  =  0.0465, 

Pe  =  =  465  for  d0  =  10"10  Pe  =  46.5  for  d0  =  10‘9  and  Pe  =  4.65 

CLq  ^  * 

2 

for  d0  =  10-8  FF  Figures  5.39  and  5.40  show  the  tracer  concentration  fields  at  four 

2  2 

times  for  the  simulations  with  d0  =  10-1°  FL.  and  d0  —  10-8  FL-  The  plots  were 
generated  directly  using  discrete  SPH  data  with  the  grey-scale  representing  the  value 
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of  concentration.  The  tracer  plume  front  becomes  obscure  and  the  concentration  field 
becomes  more  uniform  as  do  increases. 

Figure  .5.41  is  a  plot  of  M2l  as  a  function  of  time.  Values  of  M2l  were  computed 
by  subtracting  the  centered  second  moment  at  time  t  =  0  from  that  at  time  t.  Figure 
5.41  shows  that,  at  any  given  time,  tracer  spreading  about  its  center  increases  as 
d0  decreases.  Values  of  hydrodynamic  dispersion  coefficient  in  the  L  direction  DL 
were  calculated  using  Equation  5.6  and  are  shown  in  Figure  5.42.  An  asymptotic 

o 

dispersion  coefficient  was  found  for  the  simulation  with  d0  =  10-8  — .  However. 

the  curve  corresponding  to  do  =  10~10  differs  little  from  the  one  for  pure  tracer 

2 

convection.  The  curve  corresponding  to  d0  =  10"9  TJ}~  lies  between  the  curves  for 
d0  =  10"10  ^  and  d0  =  10-8  as  expected. 

Tracer  hydrodynamic  dispersion  in  two-dimensional  spatially  periodic  porous  me¬ 
dia  can  be  fundamentally  different  from  pure  tracer  convection  when  F  is  aligned 
with  a  line  of  media  symmetry.  An  asymptotic  Fickian  behavior  exists  for  tracer 
hvdrodvnamic  dispersion  under  this  condition.  However,  Fickian  behavior  appears 
some  time  after  the  tracer  is  introduced  into  the  flow.  A  characteristic  time  tc  defined 
as  (Equation  2.139), 


is  needed  for  tracer  to  sample  interstitial  space  of  the  unit  cells.  The  mechanism 

of  molecular  diffusion,  which  provides  the  possibility  for  tracer  to  sample  all  pore 

space,  is  absent  in  pure  tracer  convection.  The  value  of  tc  is  a  direct  function  of  the 

2 

coefficient  of  molecular  diffusion.  For  the  simulation  with  do  =  10  8  ,  tc  =  25  s 

if  using  Cd  =  -R  =  0.5  mm,  which  is  approximately  the  time  when  Di  becomes 

2 

essentially  constant  in  Figure  5.42.  However,  tc  —  2500  s  for  do  =  10  10  y|-  and 

2 

tc  =  250  s  for  d0  -  10~9  These  simulations  were  not  run  long  enough  to  reach  tc. 


Figure  5.39:  Tracer  concentration  fields  for  dispersion  through  a  square  array  of  circular  cylinders  (r/0  =  10  10  -£-)  at  elapsed 
times  of  (a)  Os,  (b)  33.166s,  (c.)  65.942s,  (d)  94.038s. 


Figure  5.41:  Centered  second  moment  in  the  L  direction  of  tracer  distribution  versus 
time  for  tracer  transport  through  a  square  arrav  of  circular  cylinders. 


t  (s) 


Figure  5.42:  Dispersion  coefficient  in  the  L  direction  of  tracer  distribution  versus  time 
for  tracer  transport  through  a  square  array  of  circular  cylinders. 
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5.8  Summary 

Pure  tracer  convection  and  hydrodynamic  dispersion  models  were  developed  using 
Smoothed  Particle  Hydrodynamics  (SPH).  The  dispersion  model  was  used  to  simulate 
the  Taylor  dispersion  problem  and  the  results  were  in  close  agreement  with  analytical 
solutions.  Simulations  using  the  tracer  convection  model  were  used  to  calculate  tortu¬ 
osity  and  effective  porosity  of  porous  media.  It  was  found  that  pure  tracer  convection 
through  two-dimensional  spatially  periodic  porous  media  cannot  be  described  as  an 
asymptotic  Fickian-t.vpe  process,  even  for  large  times,  if  F  is  aligned  with  a  line  of 
media  symmetry.  However,  an  asymptotic  Fickian  approximation  is  valid  for  tracer 
hydrodynamic  dispersion  through  two-dimensional  spatially  periodic  porous  media. 
If  F  is  not  aligned  with  a  line  of  media  symmetry,  Fickian-tvpe  mixing  is  possible 
for  pure  tracer  convection.  Due  to  the  time  step  constraint  for  the  explicit  time  inte¬ 
gration  method  and  limited  computer  resources,  the  SPH  model  has  limitations  for 
simulating  tracer  dispersion  through  porous  media. 
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CHAPTER  6  CONCLUSIONS  AND  RECOMMENDATIONS  FOR 

FUTURE  WORK 

6.1  Conclusions 

1.  Necessary  extensions  have  been  implemented  and  tested  which  allow  Smoothed 
Particle  Hydrodynamics  (SPH)  to  model  low  Reynolds  number  incompressible 
flow  through  porous  media: 

2.  Computed  values  of  permeability  and  cylinder  drag  force  from  simulations  of 
flow  through  spatially  periodic  porous  media  using  the  SPH  model  are  in  ac¬ 
cordance  with  Darcy’s  law  and  in  close  agreement  with  finite  element  solutions 
and  other  published  solutions  in  the  literature: 

3.  Necessary  extensions  have  been  implemented  and  tested  which  allow  SPH  to 
model  diffusion  through  porous  media; 

4.  Computed  values  of  concentration  and  effective  diffusion  coefficient  from  sim¬ 
ulations  of  diffusion  through  spatially  periodic  porous  media  using  the  SPH 
model  are  in  close  agreement  with  finite  element  solutions  and  other  published 
solutions  in  the  literature; 

5.  Nondimensional  diffusivities  d*  were  calculated  using  the  SPH  model  for  spa¬ 
tially  periodic  porous  media  with  circular  and  elliptical  cylinders  arranged  in 
square,  staggered,  and  hexagonal  arrays.  While  specific  surface  area  S0  does 
not  affect  d*,  porosity  n  and  array  type  were  found  to  be  the  most  and  least  im¬ 
portant  parameters,  respectively,  which  influence  the  values  of  d*.  Anisotropy 
of  a  medium  was  found  to  be  mainly  determined  by  the  aspect  ratio  ^  of  its 


solid  inclusions; 
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6.  Tracer  convection  and  hydrodynamic  dispersion  models  were  developed  using 
.  SPH; 

7.  Tracer  convection  simulations  were  used  to  calculate  tortuosity  and  effective 
porosity  of  porous  media.  Values  of  tortuosity  calculated  using  SPH  were  less 
than  \/2  —  1.414  as  proposed  originally  by  Carman  (1937)  for  unconsolidated 
porous  aggregates.  Tortuosity  was  found  to  be  a  function  of  media  geometry 
and  the  body  force  F  direction.  Values  of  effective  porosity  calculated  using 
SPH  were  close  to  the  nominal  porosity  values  for  the  media  type  considered  in 
this  work: 

8.  Solutions  for  Taylor  dispersion  between  infinite  plates  using  hydrodynamic  dis¬ 
persion  model  were  in  close  agreement  with  analytical  solutions; 

9.  Simulations  using  SPH  indicated  that  pure  tracer  convection  through  two- 
dimensional  spatially  periodic  porous  media  cannot  be  described  as  an  asymp¬ 
totic  Fickian-type  process,  even  for  large  times,  if  F  is  aligned  with  a  line  of 
media  symmetry.  However,  an  asymptotic  Fickian  approximation  is  valid  for 
tracer  hydrodynamic  dispersion  through  spatially  periodic  porous  media.  If  F 
is  not  aligned  with  a  line  of  media  symmetry,  Fickian-type  mixing  is  possible 
for  pure  tracer  convection; 

10.  Due  to  the  time  step  constraint  for  the  explicit  time  integration  method  and 
limited  computer  resources,  the  SPH  model  has  limitations  for  simulating  tracer 
dispersion  through  porous  media. 

6.2  Recommendations  for  Future  Work 

1.  An  advanced  parallel  algorithm  is  needed  for  the  SPH  model  to  improve  its 
efficiency  such  that  mass  transport  through  porous  media  can  be  simulated; 

2.  A  systematic  study  of  media  dispersivity  as  a  function  of  media  properties  and 
flow  conditions  can  be  conducted  using  the  SPH  model; 
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3.  Extension  of  the  SPH  model  to  three-dimensional  problems  is  needed  to  simulate 
real  random  porous  media; 

4.  Incorporation  of  progressively  more  complicated  pin  sics,  such  as  suifaco  ten¬ 
sion.  is  needed  to  model  multiphase  transport,  fluid-solid  soiption/dcsoiption 
processes,  solute  reactions,  and  decay; 

5.  A  discrete  model  for  fluid-solid  systems  can  be  achieved  by  coupling  a  three- 
dimensional  SPH  model  with  the  discrete  element  method  (DEM),  which  models 
the  solid  phase.  DEM  will  allow  modeling  of  a  deformable  solid  matrix.  It  is 
possible  to  study  the  generation  of  pore  pressure  within  porous  media  subjected 
to  both  static  and  dynamic  loadings  using  this  coupled  SPH-DEM  model.  The 
study  may  then  lead  to  important  findings  regarding  to  liquefaction  behavior 


of  loose  saturated  sand. 
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